library(data.table)
library(lfe)
library(stargazer)
library(kableExtra)
library(magrittr)
library(Matching)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(gridExtra)
library(DIDmultiplegt)
library(texreg)
library(stringr)
library(parallel)
library(sf)

# set working directory
setwd("~/Dropbox (Personal)/Trade shocks 1900s/Replication Archive")


theme_set(theme_minimal(base_family = "sans", base_size = 10))
theme_update(legend.position = "bottom")

#%% Functions

# aggregates variables to shock level and computes exposure-robust standard
# errors following Borusyak, Hull, Jaravel 2018
borusyak_reg <- function(data, y_var, x_var, shock_var, share_var, group_vars,
    control_vars = "0", fe_vars = "0") {

    # subset data from data-frame
    var_vec <- c(y_var, x_var, shock_var, share_var, group_vars)
    if (control_vars[1] != "0") {
        var_vec <- c(var_vec, control_vars)
    }
    if (fe_vars[1] != "0") {
        var_vec <- c(var_vec, fe_vars)
    }
    var_vec <- unique(var_vec)
    df <- data[, ..var_vec]


    # if control or FE variables are specified, residualize them out
    if (control_vars[1] == "0" & fe_vars[1] == "0") {
        df$y <- df[, ..y_var]
        df$x <- df[, ..x_var]
    } else {
        rx_formula <- paste0(x_var, " ~ ",
            paste(control_vars, collapse = " + "),
            " | ", paste(fe_vars, collapse = " + "))
        rx <- felm(as.formula(rx_formula), data = df)
        df[, x := rx$residuals]

        ry_formula <- paste0(y_var, " ~ ",
            paste(control_vars, collapse = " + "),
            " | ", paste(fe_vars, collapse = " + "))
        ry <- felm(as.formula(ry_formula), data = df)
        df[, y := ry$residuals]

    }
    short_vars <- c("y", "x", share_var, shock_var, group_vars)

    df_s <- df[, ..short_vars]
    colnames(df_s)[3:4] <- c("share", "shock")
    if(length(group_vars) == 2) {
        colnames(df_s)[5:6] <- c("group", "period")
    }
    ag_vars <- colnames(df_s)[4:length(colnames(df_s))]
    df_s <- df_s[, .(x = weighted.mean(x, share),
        y = weighted.mean(y, share),
        weights = sum(share)
        ),
        by = ag_vars]

    # cluster standard errors by group if multiple periods
    if (length(group_vars) < 2) {
        r_iv <- felm(y ~ 1 | 0 | (x ~ shock),
            data = df_s, weights = df_s[, weights], keepX = T)
    } else {
        r_iv <- felm(y ~ 1 | 0 | (x ~ shock) | group,
            data = df_s, weights = df_s[, weights], keepX = T)
    }
    rm(data, df, df_s)
    return(r_iv)

}


# calculates Rotemberg weights in decomposition developed by Goldsmith-Pinkham,
# Sorkin and Swift 2020
calc_rotemberg_weights <- function(data, y_var, x_var, shock_var, share_var,
    shock_group_vars, unit_group_vars, control_vars = "0", fe_vars = "0") {
    # data should be shock-level, as in shocks-elec
    var_vec <- c(y_var, x_var, shock_group_vars, unit_group_vars,
        share_var, shock_var)
    if (control_vars[1] != "0") {
        var_vec <- c(var_vec, control_vars)
    }
    if (fe_vars[1] != "0") {
        var_vec <- c(var_vec, fe_vars)
    }
    var_vec <- unique(var_vec)
    df <- data[, ..var_vec]


    df$share_id <- apply(df[, ..shock_group_vars], 1, paste, collapse = "_")
    widen_formula <- as.formula(paste0(paste(unit_group_vars, collapse = " + "),
        " ~ share_id"))
    widen_formula
    df_shares <- dcast(df, formula = widen_formula,
        value.var = share_var, fill = 0)

    outcome_vec <- c(unit_group_vars, y_var, x_var)
    if(control_vars[1] != "0") {
        outcome_vec <- c(outcome_vec, control_vars)
    }
    if(fe_vars[1] != "0"){
        outcome_vec <- c(outcome_vec, fe_vars)
    }
    outcome_vec <- unique(outcome_vec)

    df_outcomes <- unique(df[, ..outcome_vec])


    setorderv(df_outcomes, unit_group_vars)
    setorderv(df_shares, unit_group_vars)
    if (all(df_outcomes[, 1:2] == df_shares[, 1:2]) == F) {
        print("shares and outcomes do not align")
    }
    # drop missing outcome rows
    missing_vec <- as.numeric(is.na(df_outcomes[, ..y_var]))

    df_outcomes <- df_outcomes[missing_vec == 0]
    df_shares <- df_shares[missing_vec == 0]

    xr_formula <- as.formula(paste0(x_var, " ~ ",
        paste(control_vars, collapse = " + "), " | ",
        paste(fe_vars, collapse = " + ")))
    xr_formula
    X_resid <- felm(xr_formula,
        data = df_outcomes)$residuals
    yr_formula <- as.formula(paste0(y_var, " ~ ",
        paste(control_vars, collapse = " + "), " | ",
        paste(fe_vars, collapse = " + ")))
    Y_resid <- felm(yr_formula,
        data = df_outcomes)$residuals
    shock_vec <- c(shock_group_vars, "share_id", shock_var)
    df_shocks <- unique(df[, ..shock_vec])

    setorder(df_shocks, share_id)
    start_index <- length(unit_group_vars) + 1

    if(all(df_shocks$share_id ==
        colnames(df_shares)[start_index:ncol(df_shares)]) == F){
            print("shocks do not align")
        }
    G <- as.matrix(df_shocks[, ..shock_var])
    Z <- as.matrix(df_shares[, start_index:ncol(df_shares)])

    dim(G)
    dim(Z)
    dim(Y_resid)
    dim(X_resid)
    make_beta <- function(index) {
        z_k <- as.matrix(Z[, index])
        beta_k <- solve(t(z_k) %*% X_resid) *
            t(z_k) %*% Y_resid
        return(beta_k)
    }
    make_g_Z_X <- function(index) {
        z_k <- as.matrix(Z[, index])
        g_k <- G[index]
        g_Z_X <- g_k * t(z_k) %*% X_resid
        return(g_Z_X)
    }

    make_rotemberg <- function() {
        indices <- 1:nrow(G)
        bs <- sapply(indices, function(x) {make_beta(x)})
        gzks <- sapply(indices, function(x) {make_g_Z_X(x)})
        sum_gzks <- gzks
        as <- gzks / sum(gzks)
        results <- data.table(share_id = colnames(Z), beta = bs, r_weight = as)
        return(results)
    }
    df_results <- make_rotemberg()

    df_results <- merge(df_results, df_shocks, by = "share_id")
    return(df_results)
}

# extract first-stage F-stat for exposure-robust regressions
get_f_stat <- function(model) {
    f_stat <- model$stage1$iv1fstat$x[["F"]]
    f_stat <- round(f_stat, 1)
    return(f_stat)
}

#%% loading data
# main datasets for economic and political regressions:
# df_s gives constituency-by-year data, df, gives these merged with industry
# shares necessary for exposure-robust standard errors and Rotemberg weights
df <- fread("Data/_Final datasets/share_level_data.csv")

df_s <- fread("Data/_Final datasets/constituency_level_data.csv")


# news data (just contains scraped term frequencies by newspaper years)
df_news <- fread("Data/_Final datasets/newspaper_data.csv")
df_news[, term := str_replace_all(term, "\"[\"]*", "\"")]
# constituency-level data grouped to cities for cases where newspapers fall
# into multiple constituencies
df_group <- fread("Data/_Final datasets/share_level_data_grouped.csv")
df_group_s <- fread("Data/_Final datasets/constituency_level_data_grouped.csv")

# manifesto data
df_manifestos <- fread("Data/_Final datasets/manifesto_data.csv")


# functions to use with newspaper data


# aggregate observations for multiple terms by newspaper-year and standardize
aggregate_news <- function(news_df) {
    news_df <- news_df[, .(per_issue = mean(per_issue, na.rm = T)),
        by = .(constituency_group, year, title)]
    per_issue_mean <- mean(news_df$per_issue, na.rm = T)
    per_issue_sd <- sd(news_df$per_issue, na.rm = T)
    news_df[, per_issue := (per_issue - per_issue_mean) / per_issue_sd]
    news_df <- news_df[!is.na(constituency_group) & !is.na(per_issue)]
    return(news_df)
}

# regress standardized frequency on shock variable
shock_term_s <- function(news_df) {
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group_s, by = c("constituency_group", "year"),
        all.x = T, all.y = T)
    r_temp1 <- felm(per_issue ~ shock_diff_w
        | title + year | 0 |
        county_name, data = news_df)
    r_temp2 <- felm(per_issue ~ shock_diff_w + const_frac_secondary : as.factor(year)
        | title + year | 0 |
        county_name, data = news_df)
    return(list(r_temp1, r_temp2))
    rm(news_df)
}

# as above, but setting the base year for the shock variable to 1900
shock_term00s <- function(news_df) {
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group_s, by = c("constituency_group", "year"),
        all.x = T, all.y = T)
    r_temp1 <- felm(per_issue ~ shock_diff00w
        | title + year | 0 |
        county_name, data = news_df)
    r_temp2 <- felm(per_issue ~ shock_diff00w + const_frac_secondary : as.factor(year)
        | title + year | 0 |
        county_name, data = news_df)
    return(list(r_temp1, r_temp2))
    rm(news_df)
}


# regress standardized frequency on shock variable, with exposure-robust
# SEs

shock_term <- function(news_df) {
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group, by = c("constituency_group", "year"),
        all.x = T, all.y = T, allow.cartesian = T)
    r_temp1 <- borusyak_reg(data = news_df[!is.na(shock_diff_w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"))
    r_temp2 <- borusyak_reg(data = news_df[!is.na(shock_diff_w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"),
        control_vars = sec_controls_s)
    return(list(r_temp1, r_temp2))
    rm(news_df)
}
shock_term00 <- function(news_df) {
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group, by = c("constituency_group", "year"),
        all.x = T, all.y = T, allow.cartesian = T)
    r_temp1 <- borusyak_reg(data = news_df[!is.na(shock_diff00w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff00w",
        shock_var = "diff_percap00w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"))
    r_temp2 <- borusyak_reg(data = news_df[!is.na(shock_diff00w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff00w",
        shock_var = "diff_percap00w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"),
        control_vars = sec_controls_s)
    return(list(r_temp1, r_temp2))
    rm(news_df)
}

# as above, but controlling for initial non-Irish immigrant population,
# to use in regressions of anti-immigrant sentiment
shock_term_alien_s <- function(news_df){
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group_s, by = c("constituency_group", "year"),
        all.x = T, all.y = T)
    r_temp1 <- felm(per_issue ~ shock_diff_w + const_frac_n_irish : as.factor(year)
        | title + year | 0 |
        county_name, data = news_df)
    r_temp2 <- felm(per_issue ~ shock_diff_w +
        (const_frac_secondary + const_frac_n_irish) : as.factor(year)
        | title + year | 0 |
        county_name, data = news_df)
    rm(news_df)
    return(list(r_temp1, r_temp2))

}
shock_term_alien <- function(news_df) {
    news_df <- aggregate_news(news_df)

    news_df <- merge(news_df, df_group, by = c("constituency_group", "year"),
        all.x = T, all.y = T, allow.cartesian = T)
    r_temp1 <- borusyak_reg(data = news_df[!is.na(shock_diff_w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"),
        control_vars = n_irish_controls)
    r_temp2 <- borusyak_reg(data = news_df[!is.na(shock_diff_w) & !is.na(per_issue)],
        y_var = "per_issue",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("title", "year"),
        control_vars = c(n_irish_controls, sec_controls_s))
    return(list(r_temp1, r_temp2))
    rm(news_df)
}
colnames(df_group_s)

# functions to use with manifesto data
# scale by length of speech and standardize
prep_manifesto <- function(data, term) {
    # data should be df_manifestos, subset to the appropriate term and party
    data[, term_count_scaled := term_count / ntoken]
    term_count_mean <- mean(data$term_count_scaled, na.rm = T)
    term_count_sd <- sd(data$term_count_scaled, na.rm = T)
    data[, term_count_st := (term_count_scaled - term_count_mean) / term_count_sd]
    return(data)
}

# regress standardized measure on IPW, with and without MF-by-year controls,
# for Liberal candidates
manifesto_reg_lib_s <- function(data) {
    data <- prep_manifesto(data)
    data <- merge(data, df_s, by = c("election_id"), allow.cartesian = T)
    r_temp_1 <- felm(term_count_st ~ shock_diff_w | constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count_st) & party == "L"])
    r_temp_2 <- felm(term_count_st ~ shock_diff_w +
        (const_frac_secondary) : as.factor(year) |
        constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count_st) & party == "L"])
    rm(data)
    return(list(r_temp_1, r_temp_2))
}
# same but with exposure-robust SEs
manifesto_reg_lib <- function(data) {
    # data should be df_manifestos
    data <- prep_manifesto(data)
    data <- merge(data, df, by = c("election_id"), allow.cartesian = T)

    r_temp_1 <- borusyak_reg(data = data[!is.na(term_count_st) &
        !is.na(diff_percap_w) & party == "L"],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "year"),
        control_vars = "0")
    r_temp_2 <- borusyak_reg(data = data[!is.na(term_count_st) &
        !is.na(diff_percap_w) & party == "L"],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "year"),
        control_vars = c(sec_controls))
    rm(data)
    return(list(r_temp_1, r_temp_2))
}
# equivalent for Conservative candidates, controlling for initial non-Irish
# immigrants interacted with year fixed effects
manifesto_reg_alien_s <- function(data) {
    data <- prep_manifesto(data)
    data <- merge(data, df_s, by = c("election_id"), allow.cartesian = T)
    r_temp_1 <- felm(term_count_st ~ shock_diff_w +
        const_frac_n_irish : as.factor(year)| constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count_st) & party %in% c("LU", "C")])
    r_temp_2 <- felm(term_count_st ~ shock_diff_w +
        (const_frac_secondary + const_frac_n_irish) : as.factor(year) |
        constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count_st) & party %in% c("LU", "C")])
    rm(data)
    return(list(r_temp_1, r_temp_2))
}
# again, with exposure-robust SEs
manifesto_reg_alien <- function(data) {
    # data should be df_manifestos
    data <- prep_manifesto(data)
    data <- merge(data, df, by = c("election_id"), allow.cartesian = T)

    r_temp_1 <- borusyak_reg(data = data[!is.na(term_count_st) &
        !is.na(diff_percap_w) & party %in% c("LU", "C")],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "year"),
        control_vars = n_irish_controls)
    r_temp_2 <- borusyak_reg(data = data[!is.na(term_count_st) &
        !is.na(diff_percap_w) & party %in% c("LU", "C")],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "year"),
        control_vars = c(sec_controls, n_irish_controls))
    rm(data)
    return(list(r_temp_1, r_temp_2))
}



# regressions of standardized term references on IPW for all parties,
# with party FEs
df_manifestos[, party2 := party]
df_manifestos[party == "LU", party2 := "C"]

manifesto_reg_all_s <- function(data){
    data <- prep_manifesto(data)
    data <- merge(data, df_s, by = c("election_id"), allow.cartesian = T)
    r_temp_1 <- felm(term_count_st ~ shock_diff_w | party2 + constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count)])
    r_temp_2 <- felm(term_count_st ~ shock_diff_w +
        (const_frac_secondary) : as.factor(year) |
        party2 + constituency_id + year | 0 |
        county_name, data = data[!is.na(term_count)])
    rm(data)
    return(list(r_temp_1, r_temp_2))
}

manifesto_reg_all <- function(data){
    # data should be df_manifestos
    data <- prep_manifesto(data)
    data <- merge(data, df, by = c("election_id"), allow.cartesian = T)

    r_temp_1 <- borusyak_reg(data = data[!is.na(term_count) &
        !is.na(diff_percap_w)],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "party2", "year"),
        control_vars = "0")
    r_temp_2 <- borusyak_reg(data = data[!is.na(term_count) &
        !is.na(diff_percap_w)],
        y_var = "term_count_st",
        x_var = "shock_diff_w",
        shock_var = "diff_percap_w",
        share_var = "const_frac_category",
        group_vars = c("category_uk", "year"), # entity, year
        fe_vars = c("constituency_id", "party2", "year"),
        control_vars = c(sec_controls))
    rm(data)
    return(list(r_temp_1, r_temp_2))
}


# names of variables used as controls
sec_controls <- c("sec_1886", "sec_1890", "sec_1892", "sec_1895", "sec_1900",
    "sec_1906", "sec_1910", "sec_1911")
sec_controls_s <- c("sec_1886", "sec_1892", "sec_1895", "sec_1900",
        "sec_1906", "sec_1910", "sec_1911")
n_irish_controls <- c("n_irish_1886", "n_irish_1892", "n_irish_1895",
    "n_irish_1900", "n_irish_1906", "n_irish_1910", "n_irish_1911")
union_controls <- c("unions_1900", "unions_1906", "unions_1910", "unions_1911")
years <- c(1886, 1890, 1892, 1895, 1900, 1906, 1910, 1911)
steel_controls <- paste0("steel_", years)
zinc_controls <- paste0("zinc_", years)
sugar_controls <- paste0("sugar_", years)
cotton_controls <- paste0("cotton_", years)
lace_controls <- paste0("lace_", years)
pc_controls <- paste0(rep(c("pc1", "pc2", "pc3"), each = 4),
    c("_1890", "_1900", "_1906", "_1910"))

#%% regressions using economic outcome variables:

# Table 2
r_ec_1s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w |  year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

r_ec_2s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w +
    frac_unskilled_lag + frac_vagrant_lag +
    frac_secondary_lag + hiscam_lag |  year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

r_ec_3s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag + frac_vagrant_lag +
    frac_secondary_lag + hiscam_lag + const_frac_secondary : as.factor(year)|  year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

r_ec_4s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag + frac_vagrant_lag +
    frac_secondary_lag + hiscam_lag | constituency_id + year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

r_ec_5s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w  | year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

r_ec_6s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w +
    frac_unskilled_lag + frac_secondary_lag + frac_vagrant_lag +
    hiscam_lag | year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
r_ec_7s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w +
    const_frac_secondary : as.factor(year) +
    frac_unskilled_lag + frac_secondary_lag + frac_vagrant_lag +
    hiscam_lag | year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
r_ec_8s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w +
    frac_unskilled_lag + frac_secondary_lag + frac_vagrant_lag +
    hiscam_lag | constituency_id + year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])

table_econ_1s <- stargazer(r_ec_1s, r_ec_2s, r_ec_3s, r_ec_4s, r_ec_5s, r_ec_6s,
    r_ec_7s, r_ec_8s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    float.env = "sidewaystable",
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    keep = "fd_shock",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", "", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x", "", "", "", "x")),
    omit.stat = c("ser"),
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{Stacked first difference estimates,
        at the constituency level, for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects. (2)--(4) and (6)--(8) add controls
        for lagged manufacturing employment, lagged fraction in unskilled jobs,
        lagged fraction of vagrants, and lagged average economic status;
        (3) and (7) include 1880 manufacturing
        employment interacted with year dummy variables, (4) and (8) include
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by county in parentheses.}"),
    title = c("Effects of import competition on local economies"),
    label = "table_economic_effects_cl",
    notes.align = "l")

cat(table_econ_1s, sep = '\n', file = "Tables/table_economic_effects_cl.tex")

# save full regression output
table_econ_1s_full <- stargazer(r_ec_1s, r_ec_2s, r_ec_3s, r_ec_4s, r_ec_5s, r_ec_6s,
    r_ec_7s, r_ec_8s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    #keep = "fd_shock",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", "", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x", "", "", "", "x")),
    omit.stat = c("ser"),
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{Stacked first difference estimates,
        at the constituency level, for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects. (2)--(4) and (6)--(8) add controls
        for lagged manufacturing employment, lagged fraction in unskilled jobs,
        lagged fraction of vagrants, and lagged average economic status;
        (3) and (7) include 1880 manufacturing
        employment interacted with year dummy variables, (4) and (8) include
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by county in parentheses.}"),
    title = c("Full regression output for Table 2"),
    label = "table_2",
    notes.align = "l")

cat(table_econ_1s_full, sep = '\n', file = "Tables with full output/table_2.tex")

rm(r_ec_1s, r_ec_2s, r_ec_3s, r_ec_4s, r_ec_5s, r_ec_6s, r_ec_7s, r_ec_8s)

# Table A-22, replicating these results with exposure-robust standard errors
r_ec_1 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("0"))
r_ec_2 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c( "year"),
    control_vars = c("frac_unskilled_lag", "frac_vagrant_lag", "frac_secondary_lag",
        "hiscam_lag"))
r_ec_3 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c( "year"),
    control_vars = c(sec_controls, "frac_unskilled_lag", "frac_vagrant_lag",
        "frac_secondary_lag", "hiscam_lag"))
r_ec_4 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)  ],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c( "frac_unskilled_lag", "frac_vagrant_lag",
        "frac_secondary_lag", "hiscam_lag"))
r_ec_5 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("0"))
r_ec_6 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c( "year"),
    control_vars = c("frac_unskilled_lag", "frac_vagrant_lag",
        "frac_secondary_lag", "hiscam_lag"))
r_ec_7 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c( "year"),
    control_vars = c(sec_controls, "frac_unskilled_lag", "frac_vagrant_lag",
        "frac_secondary_lag", "hiscam_lag"))
r_ec_8 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)  ],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c("frac_unskilled_lag", "frac_vagrant_lag",
        "frac_secondary_lag", "hiscam_lag"))
table_econ_1 <- stargazer(r_ec_1, r_ec_2, r_ec_3, r_ec_4, r_ec_5, r_ec_6,
    r_ec_7, r_ec_8,
    dep.var.labels.include = F,
    dep.var.caption = "",
    float.env = "sidewaystable",
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", "", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x", "", "", "", "x"),
        c("First stage F-stat", get_f_stat(r_ec_1), get_f_stat(r_ec_2),
            get_f_stat(r_ec_3), get_f_stat(r_ec_4), get_f_stat(r_ec_5),
            get_f_stat(r_ec_6), get_f_stat(r_ec_7), get_f_stat(r_ec_8))),
    keep.stat = c("n"),
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table \\ref{table_economic_effects_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022).
        Stacked first difference estimates,
        at the constituency level, aggregated to the industry level,
        for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects. (2)--(4) and (6)--(8) add controls
        for lagged manufacturing employment, lagged fraction in unskilled jobs,
        lagged fraction of vagrants, and lagged average economic status;
        (3) and (7) include 1880 manufacturing
        employment interacted with year dummy variables, (4) and (8) include
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by industry in parentheses.}"),
    title = c("Effects of import competition on local economies,
        exposure-robust standard erors"),
    label = "table_economic_effects_ex",
    notes.align = "l")
cat(table_econ_1, sep = '\n', file = "Tables/table_economic_effects_ex.tex")
rm(r_ec_1, r_ec_2, r_ec_3, r_ec_4, r_ec_5, r_ec_6, r_ec_7, r_ec_8)

# Table A-2: effects on average economic status

r_ec_alt_1s <- felm(fd_hiscam ~ fd_shock_census_w |
    year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
summary(r_ec_alt_1s)


r_ec_alt_2s <- felm(fd_hiscam ~ fd_shock_census_w + hiscam_lag + frac_secondary_lag |
    year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
summary(r_ec_alt_2s)


r_ec_alt_3s <- felm(fd_hiscam ~ fd_shock_census_w + hiscam_lag +
    frac_secondary_lag + const_frac_secondary : as.factor(year) |
    year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
summary(r_ec_alt_3s)


r_ec_alt_4s <- felm(fd_hiscam ~ fd_shock_census_w + hiscam_lag + frac_secondary_lag |
    constituency_id + year | 0 | county_name,
    data = df_s[year %in% c(1890, 1900, 1910)])
summary(r_ec_alt_4s)



table_hiscam_s <- stargazer(r_ec_alt_1s, r_ec_alt_2s, r_ec_alt_3s, r_ec_alt_4s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    keep = "fd_shock",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x")),
    omit.stat = c("ser"),
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.4\\textwidth}{Stacked first difference estimates,
        at the constituency level, for 1880--1890, 1890--1900, 1900-1910.
        Dependent variable is change in average economic status.
        All models include year fixed effects. (2)--(4) add controls
        for lagged manufacturing employment and
        lagged average economic status;
        (3) includes 1880 manufacturing employment interacted with year
        dummy variables, (4) includes
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by county in parentheses.}"),
    title = c("Effects of import competition on average economic status"),
    label = "table_hiscam_cl",
    notes.align = "l")

cat(table_hiscam_s, sep = '\n',
    file = "Tables/table_hiscam_cl.tex")

table_hiscam_s_full <- stargazer(r_ec_alt_1s, r_ec_alt_2s, r_ec_alt_3s, r_ec_alt_4s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    #keep = "fd_shock",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x")),
    omit.stat = c("ser"),
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.4\\textwidth}{Stacked first difference estimates,
        at the constituency level, for 1880--1890, 1890--1900, 1900-1910.
        Dependent variable is change in average economic status.
        All models include year fixed effects. (2)--(4) add controls
        for lagged manufacturing employment and
        lagged average economic status;
        (3) includes 1880 manufacturing employment interacted with year
        dummy variables, (4) includes
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by county in parentheses.}"),
    title = c("Full regression output for Table A-2"),
    label = "table_a2",
    notes.align = "l")

cat(table_hiscam_s_full, sep = '\n',
    file = "Tables with full output/table_a2.tex")


rm(r_ec_alt_1s, r_ec_alt_2s, r_ec_alt_3s, r_ec_alt_4s)

# Table A-27: effects on average economic status, with exposure-robust SEs


r_ec_alt_1 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_hiscam",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("0"))
summary(r_ec_alt_1)

r_ec_alt_2 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_hiscam",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("hiscam_lag", "frac_secondary_lag"))
summary(r_ec_alt_2)

r_ec_alt_3 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_hiscam",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("hiscam_lag", "frac_secondary_lag",
        sec_controls))
summary(r_ec_alt_3)

r_ec_alt_4 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910) ],
    y_var = "fd_hiscam",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c("hiscam_lag", "frac_secondary_lag"))
summary(r_ec_alt_4)

table_hiscam <- stargazer(r_ec_alt_1, r_ec_alt_2, r_ec_alt_3, r_ec_alt_4,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    #keep = "fd_shock",
    omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    add.lines = list(c("Controls", "", "x", "x", "x"),
        c("Initial Mf x year", "", "", "x", ""),
        c("Constituency trends", "", "", "", "x"),
        c("First stage F-stat", get_f_stat(r_ec_alt_1), get_f_stat(r_ec_alt_2),
            get_f_stat(r_ec_alt_3), get_f_stat(r_ec_alt_4))),
    #omit.stat = c("ser"),
    keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.4\\textwidth}{
        This table replicates the results of Table \\ref{table_hiscam_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Stacked first difference estimates,
        at the constituency level, aggregated to the industry level,
        for 1880--1890, 1890--1900, 1900-1910.
        Dependent variable is change in average economic status.
        All models include year fixed effects. (2)--(4) add controls
        for lagged manufacturing employment and
        lagged average economic status;
        (3) includes 1880 manufacturing employment interacted with year
        dummy variables, (4) includes
        constituency fixed effects, which adjust for constituency-specific time
        trends.
        Standard errors clustered by industry in parentheses.}"),
    title = c("Effects of import competition on average economic status,
        exposure-robust standard errors"),
    label = "table_hiscam_ex",
    notes.align = "l")

cat(table_hiscam, sep = '\n',
    file = "Tables/table_hiscam_ex.tex")
rm(r_ec_alt_1, r_ec_alt_2, r_ec_alt_3, r_ec_alt_4)

# Table A-4: robustness checks for economic regressions
r_rob1s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + const_frac_steel : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob2s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + const_frac_zinc : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob3s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + const_frac_sugar : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob4s <- felm(fd_ln_frac_vagrant ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + (PC1 + PC2 + PC3) : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob5s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + (const_frac_steel) : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob6s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + (const_frac_zinc) : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob7s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + (const_frac_sugar) : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
r_rob8s <- felm(fd_ln_frac_unskilled ~ fd_shock_census_w + frac_unskilled_lag +
    frac_secondary_lag + frac_vagrant_lag + (PC1 + PC2 + PC3) : as.factor(year) |
    year | 0 | county_name, df_s[year %in% c(1890, 1900, 1910)])
table_econ_robust1s <- stargazer(r_rob1s, r_rob2s, r_rob3s, r_rob4s, r_rob5s, r_rob6s,
    r_rob7s, r_rob8s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    float.env = "sidewaystable",
    keep = "fd_shock",
    #omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", "x", "", "", ""),
        c("Initial zinc x year", "", "x", "", "", "", "x", "", ""),
        c("Initial sugar x year", "", "", "x", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "x", "", "", "", "x")),
    omit.stat = c("ser"),
    #keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        Stacked first difference estimates,
        at the constituency level,
        for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects, and controls for
        lagged share unskilled, lagged manufacturing employment
         and lagged fraction of vagrants;
        (1) and (5) include the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) and (6) do the same for
        employment in sheet zinc, (3) and (7) the same for sugar. (4) and (8) add
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by county in parentheses.}"),
    title = c("Robustness checks for economic variables"),
    label = "table_econ_robust_cl",
    notes.align = "l")

cat(table_econ_robust1s, sep = '\n',
    file = "Tables/table_econ_robust_cl.tex")

table_econ_robust1s_full <- stargazer(r_rob1s, r_rob2s, r_rob3s, r_rob4s, r_rob5s, r_rob6s,
    r_rob7s, r_rob8s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    #keep = "fd_shock",
    #omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", "x", "", "", ""),
        c("Initial zinc x year", "", "x", "", "", "", "x", "", ""),
        c("Initial sugar x year", "", "", "x", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "x", "", "", "", "x")),
    omit.stat = c("ser"),
    #keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        Stacked first difference estimates,
        at the constituency level,
        for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects, and controls for
        lagged share unskilled, lagged manufacturing employment
         and lagged fraction of vagrants;
        (1) and (5) include the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) and (6) do the same for
        employment in sheet zinc, (3) and (7) the same for sugar. (4) and (8) add
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by county in parentheses.}"),
    title = c("Full regression output for Table A-4"),
    label = "table_a4",
    notes.align = "l")

cat(table_econ_robust1s_full, sep = '\n',
    file = "Tables with full output/table_a4.tex")


# Table A-28: replicating robustness checks for economic variables using
# exposure-robust standard errors
r_rob1 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", steel_controls ))
r_rob2 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", zinc_controls ))
r_rob3 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", sugar_controls))
r_rob4 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", pc_controls ))
r_rob5 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", steel_controls ))
r_rob6 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", zinc_controls ))
r_rob7 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", sugar_controls ))
r_rob8 <- borusyak_reg(data = df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_unskilled",
    x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("frac_unskilled_lag", "frac_secondary_lag",
        "frac_vagrant_lag", pc_controls ))

table_econ_robust1 <- stargazer(r_rob1, r_rob2, r_rob3, r_rob4, r_rob5, r_rob6,
    r_rob7, r_rob8,
    dep.var.labels.include = F,
    dep.var.caption = "",
    float.env = "sidewaystable",
    #keep = "fd_shock",
    omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$"),
    column.labels = c("$\\Delta$ ln \\% vagrants",
        "$\\Delta$ ln \\% unskilled jobs"),
    column.separate = c(4, 4),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", "x", "", "", ""),
        c("Initial zinc x year", "", "x", "", "", "", "x", "", ""),
        c("Initial sugar x year", "", "", "x", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "x", "", "", "", "x"),
        c("First stage F-stat", get_f_stat(r_rob1), get_f_stat(r_rob2),
            get_f_stat(r_rob3), get_f_stat(r_rob4), get_f_stat(r_rob5),
            get_f_stat(r_rob6), get_f_stat(r_rob7), get_f_stat(r_rob8))),
    #omit.stat = c("ser"),
    keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        This Table replicates the results in Table \\ref{table_econ_robust_cl}
        using the aggregation and standard error calculation method recommended
        by Borusyak, Hull, and Jaravel (2022).
        Stacked first difference estimates,
        at the constituency level, aggregated to the industry level,
        for 1880--1890, 1890--1900, 1900-1910.
        All models include year fixed effects, and controls for
        lagged share unskilled, lagged manufacturing employment
         and lagged fraction of vagrants;
        (1) and (5) include the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) and (6) do the same for
        employment in sheet zinc, (3) and (7) the same for sugar. (4) and (8) add
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by industry in parentheses.}"),
    title = c("Robustness checks for economic variables,
        exposure-robust standard errors"),
    label = "table_econ_robust_ex",
    notes.align = "l")

cat(table_econ_robust1, sep = '\n',
    file = "Tables/table_econ_robust_ex.tex")
rm(r_rob1, r_rob2, r_rob3, r_rob4, r_rob5, r_rob6, r_rob7, r_rob8)

# Table A-3: Rotemberg weights for economic regressions
vagrant_rotemberg1 <- calc_rotemberg_weights(df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant", x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w", share_var = "const_frac_category",
    shock_group_vars = c("category_uk", "year"),
    unit_group_vars = c("constituency_id", "year"),
    control_vars = "0",
    fe_vars = "year")
vagrant_rotemberg1
setorder(vagrant_rotemberg1, -r_weight)
vagrant_rotemberg1 <- vagrant_rotemberg1[1:20,
    .(category_uk, year, r_weight)]
vagrant_rotemberg1


vagrant_rotemberg2 <- calc_rotemberg_weights(df[year %in% c(1890, 1900, 1910)],
    y_var = "fd_ln_frac_vagrant", x_var = "fd_shock_census_w",
    shock_var = "fd_percap_census_w", share_var = "const_frac_category",
    shock_group_vars = c("category_uk", "year"),
    unit_group_vars = c("constituency_id", "year"),
    control_vars = c("frac_unskilled_lag", "frac_vagrant_lag", "frac_secondary_lag",
        sec_controls),
    fe_vars = "year")
setorder(vagrant_rotemberg2, -r_weight)
vagrant_rotemberg2 <- vagrant_rotemberg2[1:20,
    .(category_uk, year, r_weight)]
vagrant_rotemberg2

vagrant_table <- cbind(vagrant_rotemberg1, vagrant_rotemberg2)
colnames(vagrant_table) <- rep(c("Industry", "Year", "Weight"), 2)
vagrant_table <- kable(vagrant_table, "latex", booktabs = T,
    digits = 3, caption = "Rotemberg weights for economic regressions") %>%
    add_header_above(header = c("No controls" = 3, "Controls and Mf x year" = 3)) %>%
    kable_styling(position  = "center")
vagrant_table
cat(vagrant_table, sep = "/n", file = "Tables/table_econ_rotemberg.tex")

#%% voting regressions

# Table 3: main effects on voting

# matched panel design:
# divide sample by shock from 1900 to 1910, then match constituencies on
# 1885, 1892 and 1900 Conservative voting
df_s[, shock_1910 := shock_diff00w[year == 1910], by = .(constituency_id)]

shock_1910_cutoff <- median(df_s[year == 1910, shock_1910], na.rm = T)

df_s[, treat_1910 := 1 * (shock_1910 > shock_1910_cutoff)]


df_to_match <- df_s[year %in% c(1885, 1892,  1900, 1911),
    .(constituency_id, year, cons_share, const_frac_secondary, treat_1910)]

df_to_match <- dcast(df_to_match,
    constituency_id + const_frac_secondary + treat_1910 ~ year,
    value.var = "cons_share")

df_to_match <- na.omit(df_to_match)




matches <- Match(Y = df_to_match$`1911`, Tr = df_to_match$treat_1910,
    X = df_to_match[, .(`1885`, `1892`, `1900`)],
    caliper = c(0.5, 1, 0.1))
summary(matches)

df_match <- rbind(df_to_match[matches$index.treated],
    df_to_match[matches$index.control])
df_match$weights <- rep(matches$weights, 2)

df_match <- df_match[, .(weight = sum(weights)),
    by = .(constituency_id, treat_1910)]

df_match <- merge(df_match, df_s, by = c("constituency_id", "treat_1910"),
                  all.x = T)

match_const <- unique(df_match[, .(constituency_id, weight)])
match_consts <- data.table(
    constituency_id = rep(match_const$constituency_id,
                          times = match_const$weight))

df_m <- merge(match_consts, df, by = "constituency_id", all.x = T,
              allow.cartesian = T)

# Table 3 regressions
rv_s1s <- felm(cons_share ~ shock_diff_w | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0])
rv_s2s <- felm(cons_share ~ shock_diff_w + const_frac_secondary : as.factor(year)
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0])
rv_s3s <- felm(cons_share ~ shock_diff_w | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year <= 1900])
rv_s4s <- felm(cons_share ~ shock_diff_w +
    const_frac_secondary : as.factor(year) |
    constituency_id + year | 0 | county_name,
    data = df_s[unopposed == 0 & year <= 1900])
rv_s5s <- felm(cons_share ~ shock_diff00w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900])
rv_s6s <- felm(cons_share ~ shock_diff00w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])
rv_s7s <- felm(cons_share ~ shock_diff00w |
    constituency_id + year | 0 | county_name,
    data = df_match[year >= 1900],
    weights = df_match[year >= 1900, weight])
rv_s8s <- felm(cons_share ~ shock_diff00w + const_frac_secondary : as.factor(year) |
    constituency_id + year | 0 | county_name,
    data = df_match[year >= 1900],
    weights = df_match[year >= 1900, weight])
table_votes_compact_s <- stargazer(rv_s1s, rv_s2s, rv_s3s, rv_s4s,
    rv_s5s, rv_s6s, rv_s7s, rv_s8s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "All", "All", "1885--1900", "1885--1900",
        "1900--1910", "1900-1910", "1900--1910", "1900--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("Matched panel", "", "", "", "", "", "",  "x", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2) and (4),
        (6), and (8) add manufacturing employment in 1880 interacted with
        election dummies. (7) and (8) use a panel matched on Conservative vote
        share in 1885, 1892, and 1900.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting"),
    label = "table_votes_compact_cl",
    notes.align = "l")
cat(table_votes_compact_s, sep = '\n',
    file = "Tables/table_votes_compact_cl.tex")

table_votes_compact_s_full <- stargazer(rv_s1s, rv_s2s, rv_s3s, rv_s4s,
    rv_s5s, rv_s6s, rv_s7s, rv_s8s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$", #rep(NA, 8),
        "$\\Delta\\textrm{IPW}_{1900}$"),
    order = c("shock_diff_w", "shock_diff00w"),
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "All", "All", "1885--1900", "1885--1900",
        "1900--1910", "1900-1910", "1900--1910", "1900--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("Matched panel", "", "", "", "", "", "",  "x", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2) and (4),
        (6), and (8) add manufacturing employment in 1880 interacted with
        election dummies. (7) and (8) use a panel matched on Conservative vote
        share in 1885, 1892, and 1900.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table 3"),
    label = "table_votes_compact_cl",
    notes.align = "l")
cat(table_votes_compact_s_full, sep = '\n',
    file = "Tables with full output/table_3.tex")



rm(rv_s1s, rv_s2s, rv_s3s, rv_s4s, rv_s5s, rv_s6s, rv_s7s, rv_s8s)

# Table A-23
rv_s1 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_s2 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
rv_s3 <- borusyak_reg(data = df[unopposed == 0 & year <= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_s4 <- borusyak_reg(data = df[unopposed == 0 & year <= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
rv_s5 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_s6 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
rv_s7 <- borusyak_reg(data = df_m[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = "0")
rv_s8 <- borusyak_reg(data = df_m[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
models <- list(rv_s1, rv_s2, rv_s3, rv_s4,
    rv_s5, rv_s6, rv_s7, rv_s8)
f_stats <- sapply(models, get_f_stat)


table_votes_compact <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}$"),
    omit = "Constant",
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "All", "All", "1885--1900", "1885--1900",
        "1900--1910", "1900-1910", "1900--1910", "1900--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("Matched panel", "", "", "", "", "", "",  "x", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_compact_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022)
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2) and (4),
        (6), and (8) add manufacturing employment in 1880 interacted with
        election dummies. (7) and (8) use a panel matched on Conservative vote
        share in 1885, 1892, and 1900.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting,
        exposure-robust standard errors"),
    label = "table_votes_compact_ex",
    notes.align = "l")
cat(table_votes_compact, sep = '\n',
    file = "Tables/table_votes_compact_ex.tex")
rm(models)

# Table A-5: breaking down the effect by parties
rv_1_1s <- felm(cons_share ~ shock_diff_w | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0])
rv_1_2s <- felm(cons_share ~ shock_diff_w | constituency_id + year | 0 |
    county_name,
    data = df_s[year <= 1900 & unopposed == 0])
rv_1_3s <- felm(cons_share ~ shock_diff00w | constituency_id + year | 0 |
    county_name,
    data = df_s[year >= 1900 & unopposed == 0])
rv_2s <- felm(lab_share ~ shock_diff_w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0])
rv_3_1s <- felm(lib_share ~ shock_diff_w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0])
rv_3_2s <- felm(lib_share ~ shock_diff_w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year <= 1900])
rv_3_3s <- felm(lib_share ~ shock_diff00w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900])
rv_3_4s <- felm(lib_share ~ shock_diff00w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & lab_present == 0])
table_votes_s <- stargazer(rv_1_1s, rv_1_2s, rv_1_3s, rv_2s, rv_3_1s, rv_3_2s,
    rv_3_3s, rv_3_4s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    column.labels = c("Conservative", "Labour", "Liberal vote share"),
    column.separate = c(3, 1, 4),
    add.lines = list(c("Years", "All", "1885--1900", "1900-1910", "All", "All",
        "1885-1900", "1900--1910", "1900-1910"),
        c("Excluding Labour", "", "", "", "", "", "", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for specified party.
        All models include constituency and election fixed effects.
        Model 8 excludes elections contested by Labour.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting for different parties"),
    label = "table_votes_cl",
    notes.align = "l")
cat(table_votes_s, sep = '\n', file = "Tables/table_votes_cl.tex")
rm(rv_1_1s, rv_1_2s, rv_1_3s, rv_2s, rv_3_1s, rv_3_2s, rv_3_3s, rv_3_4s)

# Table A-29 , with exposure-robust standard errors
rv_1_1 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_1_2 <- borusyak_reg(data = df[year <= 1900 & unopposed == 0],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_1_3 <- borusyak_reg(data = df[year >= 1900 & unopposed == 0 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_2 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "lab_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_3_1 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "lib_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_3_2 <- borusyak_reg(data = df[unopposed == 0 & year <= 1900 ],
    y_var = "lib_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_3_3 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900],
    y_var = "lib_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_3_4 <- borusyak_reg(
    data = df[unopposed == 0 & year >= 1900 & lab_present == 0],
    y_var = "lib_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
table_votes <- stargazer(rv_1_1, rv_1_2, rv_1_3, rv_2, rv_3_1, rv_3_2, rv_3_3,
    rv_3_4,
    covariate.labels = c("$\\Delta$IPW"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = "Constant",
    column.labels = c("Conservative", "Labour", "Liberal vote share"),
    column.separate = c(3, 1, 4),
    add.lines = list(c("Years", "All", "1885--1900", "1900-1910", "All", "All",
        "1885-1900", "1900--1910", "1900-1910"),
        c("Excluding Labour", "", "", "", "", "", "", "", "x"),
        c("First stage F-stat", get_f_stat(rv_1_1), get_f_stat(rv_1_2),
        get_f_stat(rv_1_3), get_f_stat(rv_2), get_f_stat(rv_3_1),
        get_f_stat(rv_3_2), get_f_stat(rv_3_3), get_f_stat(rv_3_4))),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Constituency-level variables aggregated to the shock level,
        exposure-robust standard errors clustered by industry in parentheses,
        all models include constituency and election fixed effects.
        Model 8 excludes elections contested by Labour.}",
    title = c("Effects of import competition on voting for different parties,
    exposure-robust
        standard errors"),
    label = "table_votes_ex",
    notes.align = "l")
cat(table_votes, sep = '\n', file = "Tables/table_votes_ex.tex")
rm(rv_1_1, rv_1_2, rv_1_3, rv_2, rv_3_1, rv_3_2, rv_3_3, rv_3_4)

# Table A-6: effects on combined Liberal and Labour vote
rv_r21s <- felm(lib_lab_share ~ shock_diff_w | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year <= 1900])
rv_r22s <- felm(lib_lab_share ~ shock_diff_w +
    const_frac_secondary : as.factor(year) |
    constituency_id + year | 0 | county_name,
    data = df_s[unopposed == 0 & year <= 1900])
rv_r23s <- felm(lib_lab_share ~ shock_diff00w  | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900])
rv_r24s <- felm(lib_lab_share ~ shock_diff00w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])
table_votes_robust_l_s <- stargazer(rv_r21s, rv_r22s, rv_r23s, rv_r24s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "1885--1900", "1885--1900", "1900--1910",
            "1900-1910"),
        c("Initial MF x election", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is
        combined share of the vote for the Liberal and Labour parties.
        All models include constituency and election fixed effects, (2) and (4)
        add the manufacturing employment in 1880 interacted with election
        dummies. Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting for combined
        Liberals and Labour"),
    label = "table_votes_robust_l_cl",
    notes.align = "l")
cat(table_votes_robust_l_s, sep = '\n', file = "Tables/table_votes_robust_l_cl.tex")

table_votes_robust_l_s_full <- stargazer(rv_r21s, rv_r22s, rv_r23s, rv_r24s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    order = c("shock_diff_w", "shock_diff00w"),
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "1885--1900", "1885--1900", "1900--1910",
            "1900-1910"),
        c("Initial MF x election", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.3\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is
        combined share of the vote for the Liberal and Labour parties.
        All models include constituency and election fixed effects, (2) and (4)
        add the manufacturing employment in 1880 interacted with election
        dummies. Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-6"),
    label = "table_a6",
    notes.align = "l")
cat(table_votes_robust_l_s_full, sep = '\n',
    file = "Tables with full output/table_a6.tex")


# Table A-30
rv_r21 <- borusyak_reg(data = df[unopposed == 0 & year <= 1900 ],
    y_var = "lib_lab_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_r22 <- borusyak_reg(data = df[unopposed == 0 & year <= 1900 ],
    y_var = "lib_lab_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
rv_r23 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "lib_lab_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
rv_r24 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "lib_lab_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
    rm(rv_r21s, rv_r22s, rv_r23s, rv_r24s)

table_votes_robust_l <- stargazer(rv_r21, rv_r22, rv_r23, rv_r24,
    covariate.labels = c("$\\Delta\\textrm{IPW}$"),
    dep.var.labels.include = F,
    #keep = "shock",
    omit = "Constant",
    dep.var.caption = "",
    add.lines = list(c("Years", "1885--1900", "1885--1900", "1900--1910", "1900-1910"),
        c("Initial MF x election", "", "x", "", "x"),
        c("First stage F-stat", get_f_stat(rv_r21), get_f_stat(rv_r22),
            get_f_stat(rv_r23), get_f_stat(rv_r24))),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        This table replicates the results of Table
        \\ref{table_votes_robust_l_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Constituency-level variables aggregated up to the industry level,
        dependent variable is combined share of the vote for the Liberal and
        Labour parties.
        All models include constituency and election fixed effects, (2) and (4) add
        the manufacturing employment in 1880 interacted with election dummies.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on voting for combined
        Liberals and Labour,
        exposure-robust standard errors"),
    label = "table_votes_robust_l_ex",
    notes.align = "l")

cat(table_votes_robust_l, sep = '\n', file = "Tables/table_votes_robust_l_ex.tex")
rm(rv_r21, rv_r22, rv_r23, rv_r24)

# A-7 voting regressions in first differences
rv_fd1s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w | year | 0 | county_name,
    data = df_s[year <= 1900 & !is.na(fd_cons_elec10)])
rv_fd2s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w +
    const_frac_secondary : as.factor(year)| year | 0 | county_name,
    data = df_s[year <= 1900 & !is.na(fd_cons_elec10)])
rv_fd3s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w | year | 0 | county_name,
    data = df_s[year > 1900 & !is.na(fd_cons_elec10)])
rv_fd4s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w + const_frac_secondary :
    as.factor(year) | year | 0 | county_name,
    data = df_s[year > 1900 & !is.na(fd_cons_elec10)])
rv_fd5s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w + fd_shock_elec10_w_lead |
    year | 0 | county_name,
    data = df_s[year <= 1900 & !is.na(fd_cons_elec10)])
rv_fd6s <- felm(fd_cons_elec10 ~ fd_shock_elec10_w + fd_shock_elec10_w_lead +
    const_frac_secondary : as.factor(year)| year | 0 | county_name,
    data = df_s[year <= 1900 & !is.na(fd_cons_elec10)])
table_votes_fd_s <- stargazer(rv_fd1s, rv_fd2s, rv_fd3s, rv_fd4s, rv_fd5s,
    rv_fd6s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$",
        "$\\Delta\\textrm{IPW}_{t + 1}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years", "1885--1900", "1885--1900", "1900--1910",
        "1900-1910", "1885--1900", "1885--1900"),
        c("Initial Mf x election", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        Constituency-level stacked first-difference regressions, for waves
        1885--1892, 1892--1900, 1900-1910 (note there were two elections in 1910).
        Dependent variable is change in share of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2), (4) and
        (6) add controls for manufacturing employment in 1880 interacted with election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("First-difference effects of import competition on voting"),
    label = "table_votes_fd_cl",
    notes.align = "l")
cat(table_votes_fd_s, sep = '\n', file = "Tables/table_votes_fd_cl.tex")
table_votes_fd_s_full <- stargazer(rv_fd1s, rv_fd2s, rv_fd3s, rv_fd4s, rv_fd5s,
    rv_fd6s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_t$",
        "$\\Delta\\textrm{IPW}_{t + 1}$"),
    dep.var.labels.include = F,
    #keep = "shock",
    order = c("fd_shock_elec10_w", "fd_shock_elec10_w_lead"),
    dep.var.caption = "",
    add.lines = list(c("Years", "1885--1900", "1885--1900", "1900--1910",
        "1900-1910", "1885--1900", "1885--1900"),
        c("Initial Mf x election", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.3\\textwidth}{
        Constituency-level stacked first-difference regressions, for waves
        1885--1892, 1892--1900, 1900-1910 (note there were two elections in 1910).
        Dependent variable is change in share of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2), (4) and
        (6) add controls for manufacturing employment in 1880 interacted with election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-7"),
    label = "table_a7",
    notes.align = "l")
cat(table_votes_fd_s_full, sep = '\n', file = "Tables with full output/table_a7.tex")



rm(rv_fd1s, rv_fd2s, rv_fd3s, rv_fd4s, rv_fd5s, rv_fd6s)

# Table A-31: first-difference voting regressions with exposure-robust SEs
rv_fd1 <- borusyak_reg(data = df[year <= 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w",
    shock_var = "fd_percap_elec10_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = "0")
rv_fd2 <- borusyak_reg(data = df[ year <= 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w",
    shock_var = "fd_percap_elec10_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = sec_controls)
rv_fd3 <- borusyak_reg(data = df[year > 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w",
    shock_var = "fd_percap_elec10_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = "0")
rv_fd4 <- borusyak_reg(data = df[year > 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w",
    shock_var = "fd_percap_elec10_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = sec_controls)
rv_fd5 <- borusyak_reg(data = df[year <= 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w_lead",
    shock_var = "fd_percap_elec10_w_lead",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = "fd_shock_elec10_w")
rv_fd6 <- borusyak_reg(data = df[year <= 1900 & !is.na(fd_cons_elec10)],
    y_var = "fd_cons_elec10",
    x_var = "fd_shock_elec10_w_lead",
    shock_var = "fd_percap_elec10_w_lead",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"),
    control_vars = c("fd_shock_elec10_w", sec_controls))
table_votes_fd <- stargazer(rv_fd1, rv_fd2, rv_fd3, rv_fd4, rv_fd5,
    rv_fd6,
    covariate.labels = c("Shock"),
    dep.var.labels.include = F,
    #keep = "shock",
    omit = "Constant",
    dep.var.caption = "",
    add.lines = list(c("Shock variable", "$\\Delta\\textrm{IPW}_t$",
        "$\\Delta\\textrm{IPW}_t$", "$\\Delta\\textrm{IPW}_t$",
        "$\\Delta\\textrm{IPW}_t$", "$\\Delta\\textrm{IPW}_{t+1}$",
        "$\\Delta\\textrm{IPW}_{t+1}$"),
        c("Years", "1885--1900", "1885--1900", "1900--1910",
        "1900-1910", "1885--1900", "1885--1900"),
        c("Initial Mf x election", "", "x", "", "x", "", "x"),
        c("First stage F stat", get_f_stat(rv_fd1), get_f_stat(rv_fd2),
            get_f_stat(rv_fd3), get_f_stat(rv_fd4), get_f_stat(rv_fd5),
            get_f_stat(rv_fd6))),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        This table replicates the results of Table
        \\ref{table_votes_fd_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022).
        Stacked first-difference regressions, for waves
        1885--1892, 1892--1900, 1900-1910 (note there were two elections in 1910),
        with constituency-level variables aggregated up to the industry level.
        Dependent variable is change in share of the vote for the Conservative Party.
        All models include constituency and election fixed effects, (2), (4) and
        (6) add controls for manufacturing employment in 1880 interacted with election dummies.
        Standard errors clustered by industry in parentheses, (5) and (6) control
        for $\\Delta$ IPW for the correct period.}",
    title = c("First-difference effects of import competition on voting,
        exposure-robust standard errors"),
    label = "table_votes_fd_ex",
    notes.align = "l")
cat(table_votes_fd, sep = '\n', file = "Tables/table_votes_fd_ex.tex")
rm(rv_fd1, rv_fd2, rv_fd3, rv_fd4, rv_fd5, rv_fd6)

# Table A-8: Rotemberg weights for voting regressions
voting00_rotemberg1 <- calc_rotemberg_weights(df[year >= 1900 & unopposed == 0],
    y_var = "cons_share", x_var = "shock_diff00w",
    shock_var = "diff_percap00w", share_var = "const_frac_category",
    shock_group_vars = c("category_uk", "year"),
    unit_group_vars = c("constituency_id", "year"),
    control_vars = c("0"),
    fe_vars = c("constituency_id", "year"))
setorder(voting00_rotemberg1, -r_weight)
voting00_rotemberg1 <- voting00_rotemberg1[1:20, .(category_uk, year, r_weight)]

voting00_rotemberg2 <- calc_rotemberg_weights(df[year >= 1900 & unopposed == 0],
    y_var = "cons_share", x_var = "shock_diff00w",
    shock_var = "diff_percap00w", share_var = "const_frac_category",
    shock_group_vars = c("category_uk", "year"),
    unit_group_vars = c("constituency_id", "year"),
    control_vars = c(sec_controls),
    fe_vars = c("constituency_id", "year"))
setorder(voting00_rotemberg2, -r_weight)
voting00_rotemberg2 <- voting00_rotemberg2[1:20, .(category_uk, year, r_weight)]

voting_rotem_table <- cbind(voting00_rotemberg1, voting00_rotemberg2)
colnames(voting_rotem_table) <- rep(c("Industry", "Year", "Weight"), 2)
voting_rotem_table <- kable(voting_rotem_table, "latex", booktabs = T,
    digits = 3, caption = "Rotemberg weights for post-1900 voting regressions") %>%
    add_header_above(header = c("No controls" = 3, "Initial Mf x election" = 3)) %>%
    kable_styling(position  = "center")
voting_rotem_table
cat(voting_rotem_table, sep = "/n", file = "Tables/table_voting_rotemberg.tex")

# Table A-9: voting regressions controlling for different initial industries
r_elec_rob_1s <- felm(cons_share ~ shock_diff00w +
    as.factor(year) : (const_frac_secondary + const_frac_steel) |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900])
r_elec_rob_2s <- felm(cons_share ~ shock_diff00w +
    as.factor(year) : (const_frac_secondary + const_frac_cotton) |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900])
r_elec_rob_3s <- felm(cons_share ~ shock_diff00w +
    as.factor(year) : (const_frac_secondary + const_frac_sugar) |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900])
r_elec_rob_4s <- felm(cons_share ~ shock_diff00w +
    as.factor(year) : (const_frac_secondary + const_frac_lace) |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900])
r_elec_rob_5s <- felm(cons_share ~ shock_diff00w +
    as.factor(year) : (const_frac_secondary +  PC1 + PC2 + PC3) |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900])
table_elec_robust1s <- stargazer(r_elec_rob_1s, r_elec_rob_2s, r_elec_rob_3s,
    r_elec_rob_4s, r_elec_rob_5s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    keep = "shock",
    #omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1900}$"),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", ""),
        c("Initial cotton x year", "", "x", "", "", ""),
        c("Initial sugar x year", "", "", "x", "", ""),
        c("Initial lace x year", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "", "x")),
    omit.stat = c("ser"),
    #keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regressions,
        for 1900--1910. Dependent variable is share of the vote for
        Conservative candidates.
        All models include constituency and year fixed effects, and initial
        manufacturing by year controls.
        (1) includes the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) does the same for
        employment in sheet zinc, (3) does the same for sugar,
        (4) does the same for lace. (5) adds
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by county in parentheses.}"),
    title = c("Robustness checks for post-1900 voting regressions"),
    label = "table_elec_robust_cl",
    notes.align = "l")

cat(table_elec_robust1s, sep = '\n',
    file = "Tables/table_elec_robust_cl.tex")

table_elec_robust1s_full <- stargazer(r_elec_rob_1s, r_elec_rob_2s, r_elec_rob_3s,
    r_elec_rob_4s, r_elec_rob_5s,
    dep.var.labels.include = F,
    dep.var.caption = "",
    font.size = "scriptsize",
    #float.env = "sidewaystable",
    #keep = "shock",
    #omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1900}$"),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", ""),
        c("Initial cotton x year", "", "x", "", "", ""),
        c("Initial sugar x year", "", "", "x", "", ""),
        c("Initial lace x year", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "", "x")),
    omit.stat = c("ser"),
    #keep.stat = "n",
    #column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.3\\textwidth}{
        Constituency-level fixed effects regressions,
        for 1900--1910. Dependent variable is share of the vote for
        Conservative candidates.
        All models include constituency and year fixed effects, and initial
        manufacturing by year controls.
        (1) includes the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) does the same for
        employment in sheet zinc, (3) does the same for sugar,
        (4) does the same for lace. (5) adds
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by county in parentheses.}"),
    title = c("Full regression output for Table A-9"),
    label = "table_elec_robust_cl",
    notes.align = "l")

cat(table_elec_robust1s_full, sep = '\n',
    file = "Tables with full output/table_a9.tex")



rm(r_elec_rob_1s, r_elec_rob_2s, r_elec_rob_3s, r_elec_rob_4s, r_elec_rob_5s)

# Table A-32: redone with exposure-robust SEs
r_elec_rob_1 <- borusyak_reg(data = df[year >= 1900 & unopposed == F],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(steel_controls, sec_controls))
r_elec_rob_2 <- borusyak_reg(data = df[year >= 1900 & unopposed == F],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(cotton_controls, sec_controls))
r_elec_rob_3 <- borusyak_reg(data = df[year >= 1900 & unopposed == F],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(sugar_controls, sec_controls))
r_elec_rob_4 <- borusyak_reg(data = df[year >= 1900 & unopposed == F],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(lace_controls, sec_controls))
r_elec_rob_5 <- borusyak_reg(data = df[year >= 1900 & unopposed == F],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(pc_controls, sec_controls))
table_elec_robust1 <- stargazer(r_elec_rob_1, r_elec_rob_2, r_elec_rob_3,
    r_elec_rob_4, r_elec_rob_5,
    dep.var.labels.include = F,
    dep.var.caption = "",
    #float.env = "sidewaystable",
    #keep = "fd_shock",
    omit = "Constant",
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1900}$"),
    add.lines = list(
        c("Initial steel x year", "x", "", "", "", ""),
        c("Initial cotton x year", "", "x", "", "", ""),
        c("Initial sugar x year", "", "", "x", "", ""),
        c("Initial lace x year", "", "", "", "x", ""),
        c("Initial shares PCA x year", "", "", "", "", "x"),
        c("First stage F-stat", get_f_stat(r_elec_rob_1),
            get_f_stat(r_elec_rob_2), get_f_stat(r_elec_rob_3),
            get_f_stat(r_elec_rob_4), get_f_stat(r_elec_rob_5))),
    #omit.stat = c("ser"),
    keep.stat = "n",
    column.sep.width = "-10pt",
    notes = c("\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table \\ref{table_elec_robust_cl}
        using the aggregation and standard error calculation methods recommended
        by Borusyak, Hull, and Jaravel (2022).
        Constituency-level fixed effects regressions, aggregated to the
        industry level for exposure-robust standard errors,
        for 1900--1910. Dependent variable is share of the vote for
        Conservative candidates.
        All models include constituency and year fixed effects, and initial
        manufacturing by year controls.
        (1) includes the share of employment in 1881 in sheet iron and
        steel interacted with year fixed effects, (2) does the same for
        employment in sheet zinc, (3) does the same for sugar,
        (4) does the same for lace. (5) adds
        the first three principal components for the 1881 industry shares
        interacted with year fixed effects.
        Standard errors clustered by industry in parentheses.}"),
    title = c("Robustness checks for post-1900 voting regressions,
        exposure-robust standard errors"),
    label = "table_elec_robust_ex",
    notes.align = "l")

cat(table_elec_robust1, sep = '\n',
    file = "Tables/table_elec_robust_ex.tex")
rm(r_elec_rob_1, r_elec_rob_2, r_elec_rob_3, r_elec_rob_4, r_elec_rob_5)

# A-10: regressions of Conservative voting on IPW dropping years
rvd1s <- felm(cons_share ~ shock_diff00w
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1900])
rvd2s <- felm(cons_share ~ shock_diff00w + const_frac_secondary : as.factor(year)
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1900])

rvd3s <- felm(cons_share ~ shock_diff00w
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1906])
rvd4s <- felm(cons_share ~ shock_diff00w + const_frac_secondary : as.factor(year)
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1906])

rvd5s <- felm(cons_share ~ shock_diff00w
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1910])
rvd6s <- felm(cons_share ~ shock_diff00w + const_frac_secondary : as.factor(year)
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1910])

rvd7s <- felm(cons_share ~ shock_diff00w
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1911])
rvd8s <- felm(cons_share ~ shock_diff00w + const_frac_secondary : as.factor(year)
    | constituency_id + year | 0 |
    county_name,
    data = df_s[unopposed == 0 & year >= 1900 & year != 1911])

models <- list(rvd1s, rvd2s, rvd3s, rvd4s, rvd5s, rvd6s, rvd7s, rvd8s)
table_votes_drop_s <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Excluding", "1900", "1900", "1906", "1906", "1910J",
        "1910J", "1910D", "1910D"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Each model drops one election from the period.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting, dropping certain years"),
    label = "table_votes_drop_cl",
    notes.align = "l")
cat(table_votes_drop_s, sep = '\n',
    file = "Tables/table_votes_drop_cl.tex")

table_votes_drop_s_full <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Excluding", "1900", "1900", "1906", "1906", "1910J",
        "1910J", "1910D", "1910D"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    #column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Each model drops one election from the period.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression results for Table A-10"),
    label = "table_votes_drop_cl",
    notes.align = "l")
cat(table_votes_drop_s_full, sep = '\n',
    file = "Tables with full output/table_a10.tex")


rm(models)

# Table A-33
rvd1 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1900],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = "0")
rvd2 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1900],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
summary(rvd2)
rvd3 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1906],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = "0")
rvd4 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1906],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
summary(rvd4)
rvd5 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1910],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = "0")
rvd6 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1910],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
summary(rvd6)
rvd7 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1911],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = "0")
rvd8 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 & year != 1911],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = sec_controls)
summary(rvd8)
models <- list(rvd1, rvd2, rvd3, rvd4, rvd5, rvd6, rvd7, rvd8)
f_stats <- sapply(models, get_f_stat)

table_votes_drop <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    omit = "Constant",
    dep.var.caption = "",
    add.lines = list(c("Excluding", "1900", "1900", "1906", "1906", "1910J",
        "1910J", "1910D", "1910D"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_drop_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022).
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Each model drops one election from the period.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on voting, dropping certain years,
        exposure-robust standard errors"),
    label = "table_votes_drop_ex",
    notes.align = "l")
cat(table_votes_drop, sep = '\n',
    file = "Tables/table_votes_drop_ex.tex")

# A-11: moderating effects of unions on relationship between IPW and voting
median_unionists <- median(df_s[, webb_frac_unionists_1892], na.rm = T)
ru1s <- felm(cons_share ~ shock_diff_w / webb_frac_unionists_1892  | constituency_id + year | 0 | webb_county,
    df_s[])

ru2s <- felm(cons_share ~ shock_diff_w / webb_frac_unionists_1892 +
    as.factor(year) : const_frac_secondary | constituency_id + year | 0 | webb_county,
    df_s[])

ru3s <- felm(cons_share ~ shock_diff00w / webb_frac_unionists_1892  |
    constituency_id + year | 0 | webb_county,
    df_s[year >= 1900])

ru4s <- felm(cons_share ~ shock_diff00w / webb_frac_unionists_1892 +
    as.factor(year) : const_frac_secondary | constituency_id + year | 0 | webb_county,
    df_s[year >= 1900])

ru5s <- felm(cons_share ~ shock_diff_w  + as.factor(year) : const_frac_secondary | constituency_id + year | 0 | webb_county,
    df_s[webb_frac_unionists_1892 > median_unionists])

ru6s <- felm(cons_share ~ shock_diff_w  + as.factor(year) : const_frac_secondary | constituency_id + year | 0 | webb_county,
    df_s[webb_frac_unionists_1892 <= median_unionists])

ru7s <- felm(cons_share ~ shock_diff00w  + as.factor(year) : const_frac_secondary
    | constituency_id + year | 0 | webb_county,
    df_s[webb_frac_unionists_1892 > median_unionists
        & year >= 1900])

ru8s <- felm(cons_share ~ shock_diff00w  + as.factor(year) : const_frac_secondary
    | constituency_id + year | 0 | webb_county,
    df_s[webb_frac_unionists_1892 <= median_unionists &
        year >= 1900])

ru9s <- felm(cons_share ~ shock_diff00w + as.factor(year) : (webb_frac_unionists_1892)
    | constituency_id + year | 0 | webb_county,
    df_s[year >= 1900])

ru10s <- felm(cons_share ~ shock_diff00w + as.factor(year) : (const_frac_secondary
    + webb_frac_unionists_1892) | constituency_id + year | 0 | webb_county,
    df_s[year >= 1900])




table_union_s <- stargazer(ru1s, ru2s, ru3s, ru4s, ru5s, ru6s, ru7s,
    ru8s, ru9s, ru10s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1885} \\times \\textrm{\\% Union}$",
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{IPW}_{1900} \\times \\textrm{\\% Union}$"),
    dep.var.labels.include = F,
    omit = "as.factor",
    dep.var.caption = "",
    add.lines = list(c("Years", "All", "All", "1900--1910", "1900--1910",
    "All", "All", "1900--1910", "1900--1910", "1900--1910", "1900--1910"),
        c("Union sub-sample", "All", "All", "All", "All", "2H",
        "1H", "2H", "1H", "All", "All"),
        c("Union x election", "", "", "", "", "", "", "", "", "x", "x"),
        c("Initial MF x election", "", "x", "", "x", "x", "x", "x", "x", "", "x")
        ),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.9\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party. Data on union membership relative
        to population in 1892 at the county level is taken from Sidney and Beatrice Webb,
        \\emph{The History of Trade Unionism} (London: Longmans, Green and Co., 1896).
        Models (5) and (7) are estimated for constituencies with above-median
        unionization, (6) and (8) for constituencies with below-median unionization.
        Models (9) and (10) replicate regressions from table \\ref{table_votes_compact_cl},
        adding controls for unionization interacted with year dummy variables.
        All models include constituency and election fixed effects,
        Standard errors clustered by county in parentheses.}",
    title = c("Moderating effect of unions on effect of import competition on voting"),
    label = "table_union_cl",
    notes.align = "l")
cat(table_union_s, sep = '\n',
    file = "Tables/table_union_cl.tex")

table_union_s_full <- stargazer(ru1s, ru2s, ru3s, ru4s, ru5s, ru6s, ru7s,
    ru8s, ru9s, ru10s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1885} \\times \\textrm{\\% Union}$",
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{IPW}_{1900} \\times \\textrm{\\% Union}$"),
    order = c("shock_diff_w", "shock_diff_w:webb_frac_unionists_1892",
        "shock_diff00w", "shock_diff00w:webb_frac_unionists_1892")   ,
    dep.var.labels.include = F,
    #omit = "as.factor",
    dep.var.caption = "",
    add.lines = list(c("Years", "All", "All", "1900--1910", "1900--1910",
    "All", "All", "1900--1910", "1900--1910", "1900--1910", "1900--1910"),
        c("Union sub-sample", "All", "All", "All", "All", "2H",
        "1H", "2H", "1H", "All", "All"),
        c("Union x election", "", "", "", "", "", "", "", "", "x", "x"),
        c("Initial MF x election", "", "x", "", "x", "x", "x", "x", "x", "", "x")
        ),
    font.size = "scriptsize",
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party. Data on union membership relative
        to population in 1892 at the county level is taken from Sidney and Beatrice Webb,
        \\emph{The History of Trade Unionism} (London: Longmans, Green and Co., 1896).
        Models (5) and (7) are estimated for constituencies with above-median
        unionization, (6) and (8) for constituencies with below-median unionization.
        Models (9) and (10) replicate regressions from table \\ref{table_votes_compact_cl},
        adding controls for unionization interacted with year dummy variables.
        All models include constituency and election fixed effects,
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-11"),
    label = "table_union_cl",
    notes.align = "l")
cat(table_union_s_full, sep = '\n',
    file = "Tables with full output/table_a11.tex")


rm(ru1s, ru2s, ru3s, ru4s, ru5s, ru6s, ru7s,
        ru8s, ru9s, ru10s)

# A-34: moderating effects of unions with exposure robust SEs
ru5 <- borusyak_reg(data = df[unopposed == 0 & webb_frac_unionists_1892 > median_unionists],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(sec_controls))


ru6 <- borusyak_reg(data = df[unopposed == 0 & webb_frac_unionists_1892 <= median_unionists],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(sec_controls))

ru7 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 &
    webb_frac_unionists_1892 > median_unionists],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(sec_controls))

ru8 <- borusyak_reg(data = df[unopposed == 0 &
    webb_frac_unionists_1892 <= median_unionists & year >= 1900],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(sec_controls))

ru9 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(union_controls))

ru10 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"),
    control_vars = c(union_controls, sec_controls))

models <- list(ru5, ru6, ru7, ru8, ru9, ru10)
f_stats <- sapply(models, get_f_stat)

table_union <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}$"),
    omit = "Constant",
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Years",
    "All", "All", "1900--1910", "1900--1910", "1900--1910", "1900--1910"),
        c("Union sub-sample", "2H",
        "1H", "2H", "1H", "All", "All"),
        c("Union x election",  "", "", "", "", "x", "x"),
        c("Initial MF x election",  "x", "x", "x", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        This table replicates the results of models (5)--(10) of Table \\ref{table_union_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022) (their method does not allow us to estimate standard errors for
        variables interacted with the shock).
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party. Data on union membership relative
        to population in 1892 at the county level is taken from Sidney and Beatrice Webb,
        \\emph{The History of Trade Unionism} (London: Longmans, Green and Co., 1896).
        Models (1) and (3) are estimated for constituencies with above-median
        unionization, (2) and (4) for constituencies with below-median unionization.
        Models (5) and (6) replicate regressions from table \\ref{table_votes_compact_cl},
        adding controls for unionization interacted with year dummy variables.
        All models include constituency and election fixed effects,
        Standard errors clustered by industry in parentheses.}",
    title = c("Moderating effect of unions on effect of import competition on voting,
        exposure-robust standard errors"),
    label = "table_union_ex",
    notes.align = "l")
cat(table_union, sep = '\n',
    file = "Tables/table_union_ex.tex")

# A-12: Effects of imports on voting broken down by period
r_period1s <- felm(cons_share ~ shock_diff_w | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1885, 1895)])
r_period2s <- felm(cons_share ~ shock_diff_w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1885, 1895)])
r_period3s <- felm(cons_share ~ shock_diff00w | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1895, 1906)])
r_period4s <- felm(cons_share ~ shock_diff00w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1895, 1906)])

r_period5s <- felm(cons_share ~ shock_diff00w | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1906, 1911)])
r_period6s <- felm(cons_share ~ shock_diff00w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1906, 1911)])
r_period7s <- felm(cons_share ~ shock_diff00w | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1895, 1911)])
r_period8s <- felm(cons_share ~ shock_diff00w +
    const_frac_secondary : as.factor(year) | constituency_id + year | 0
    | county_name, df_s[unopposed == 0 & between(year, 1895, 1911)])
models <- list(r_period1s, r_period2s, r_period3s, r_period4s, r_period5s,
    r_period6s,  r_period7s, r_period8s)
table_votes_by_period_s <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    column.labels = c("1885--1895", "1895--1906", "1906--1910",
        "1895--1910"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(#c("Years", "1885--1892", "1885--1892", "1892--1900",
        #"1892--1900", "1900--1910", "1900--1910", "1892--1895", "1892--1895",
        #"1895--1910", "1895--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, subset by different groups of years.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting, by period"),
    label = "table_votes_by_period_cl",
    notes.align = "l")
cat(table_votes_by_period_s, sep = '\n',
    file = "Tables/table_votes_by_period_cl.tex")

table_votes_by_period_s_full <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    #keep = "shock",
    order = c("shock_diff_w", "shock_diff00w"),
    dep.var.caption = "",
    column.labels = c("1885--1895", "1895--1906", "1906--1910",
        "1895--1910"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(#c("Years", "1885--1892", "1885--1892", "1892--1900",
        #"1892--1900", "1900--1910", "1900--1910", "1892--1895", "1892--1895",
        #"1895--1910", "1895--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.4\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, subset by different groups of years.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-12"),
    label = "table_votes_by_period_cl",
    notes.align = "l")
cat(table_votes_by_period_s_full, sep = '\n',
    file = "Tables with full output/table_a12.tex")


rm(r_period1s, r_period2s, r_period3s, r_period4s, r_period5s,
        r_period6s,  r_period7s, r_period8s)
rm(models)

# A-35: effects on voting by period, exposure-robust SEs
r_period1 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1885, 1895) ],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    #control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_period2 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1885, 1895) ],
    y_var = "cons_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_period3 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1895, 1906) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    #control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_period4 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1895, 1906) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))


r_period5 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1906, 1911) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    #control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_period6 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1906, 1911) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))


r_period7 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1895, 1911) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    #control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_period8 <- borusyak_reg(data = df[unopposed == 0 & between(year, 1895, 1911) ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    control_vars = sec_controls,
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))


models <- list(r_period1, r_period2, r_period3, r_period4, r_period5,
    r_period6, r_period7, r_period8)
f_stats <- sapply(models, get_f_stat)
table_votes_by_period <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}$"),
    dep.var.labels.include = F,
    keep = "x",
    dep.var.caption = "",
    column.labels = c("1885--1895", "1895--1906", "1906--1910", "1895--1910"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(#c("Years", "1885--1892", "1885--1892", "1892--1900",
        #"1892--1900", "1900--1910", "1900--1910", "1892--1895", "1892--1895",
        #"1895--1910", "1895--1910"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_by_period_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022)
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party, subset by different groups of years.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on voting, by period, with
        exposure-robust standard errors"),
    label = "table_votes_by_period_ex",
    notes.align = "l")
cat(table_votes_by_period, sep = '\n',
    file = "Tables/table_votes_by_period_ex.tex")
rm(r_period1, r_period2, r_period3, r_period4, r_period5,
        r_period6, r_period7, r_period8)

# A-13: effects of IPW on voting, controlling for exports and wheat imports
r_r3_1s <- felm(cons_share ~ shock_diff00w + export_shock00w | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])

r_r3_2s <- felm(cons_share ~ shock_diff00w + export_shock00w +
        as.factor(year) : const_frac_secondary | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])

r_r3_3s <- felm(cons_share ~ shock_diff00w + wheat_us_shock | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])

r_r3_4s <- felm(cons_share ~ shock_diff00w + wheat_us_shock +
    as.factor(year) : const_frac_secondary | constituency_id + year | 0 |
        county_name,
    data = df_s[unopposed == 0 & year >= 1900])

models <- list(r_r3_1s, r_r3_2s, r_r3_3s, r_r3_4s)
table_votes_rev1_s <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{Exports per worker}_{1900}$",
        "$\\Delta\\textrm{US wheat imports per worker}_{1900}$"),
    dep.var.labels.include = F,
    keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Initial MF x election", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.45\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Models 1 and 2 in addition control for exports to Germany
        per worker, computed the same way as $\\Delta\\textrm{IPW}$, models 3 and
        4 control for US wheat imports per worker, with wheat employment calculated
        using agricultural laborers weighted by the share of county land devoted
        to wheat cultivation.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on voting, controlling for exports and wheat imports"),
    label = "table_votes_rev1_cl",
    notes.align = "l")
cat(table_votes_rev1_s, sep = '\n',
    file = "Tables/table_votes_rev1_cl.tex")

table_votes_rev1_s_full <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{Exports per worker}_{1900}$",
        "$\\Delta\\textrm{US wheat imports per worker}_{1900}$"),
    order = c("shock_diff00w", "export_shock00w", "wheat_us_shock"),
    dep.var.labels.include = F,
    #keep = "shock",
    dep.var.caption = "",
    add.lines = list(c("Initial MF x election", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.3\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Models 1 and 2 in addition control for exports to Germany
        per worker, computed the same way as $\\Delta\\textrm{IPW}$, models 3 and
        4 control for US wheat imports per worker, with wheat employment calculated
        using agricultural laborers weighted by the share of county land devoted
        to wheat cultivation.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-13"),
    label = "table_votes_rev1_cl",
    notes.align = "l")
cat(table_votes_rev1_s_full, sep = '\n',
    file = "Tables with full output/table_a13.tex")


# A-36 effects of IPW on voting controlling for export and wheat shocks
# with exposure-robust SEs
r_r3_1 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = "export_shock00w",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_r3_2 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("export_shock00w", sec_controls),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_r3_3 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("wheat_us_shock"),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_r3_4 <- borusyak_reg(data = df[unopposed == 0 & year >= 1900 ],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("wheat_us_shock", sec_controls),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

models <- list(r_r3_1, r_r3_2, r_r3_3, r_r3_4)
f_stats <- sapply(models, get_f_stat)

table_votes_rev1 <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    keep = "x",
    dep.var.caption = "",
    add.lines = list(c("Initial MF x election", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.45\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_rev1_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022)
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party, for the period 1900--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies. Models 1 and 2 in addition control for exports to Germany
        per worker, computed the same way as $\\Delta\\textrm{IPW}$, models 3 and
        4 control for US wheat imports per worker, with wheat employment calculated
        using agricultural laborers weighted by the share of county land devoted
        to wheat cultivation.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on voting, controlling for
        exports and wheat imports, with exposure-robust standard errors"),
    label = "table_votes_rev1_ex",
    notes.align = "l")
cat(table_votes_rev1, sep = '\n',
    file = "Tables/table_votes_rev1_ex.tex")

# A-14: regressions of IPW on voting, checks for pre-trends
# computing conservative vote predicted by constituency-specific trends
constituencies <- unique(df_s[!is.na(constituency_id) & !is.na(cons_share), constituency_id])


# predicted Conservative vote based on constituency trend 1885--1910
make_pred_cons <- function(my_constituency) {
    data_subset <- df_s[constituency_id == my_constituency &
        unopposed == 0]
    if (nrow(data_subset) == 0) {
        return("No data")
    } else {
        model <- lm(cons_share ~ year, data_subset)
        df_s[constituency_id == my_constituency,
            cons_share_pred := model$coef[1] + model$coef[2] * year]
    }


}

# predicted Conservative vote based on constituency trend 1885--1900
make_pred_cons_pre <- function(my_constituency) {
    data_subset <- df_s[constituency_id == my_constituency &
        unopposed == 0 & year <= 1900]
    if (nrow(data_subset) == 0) {
        return("No data")
    } else {
        model <- lm(cons_share ~ year, data_subset)
        df_s[constituency_id == my_constituency,
            cons_share_pred_pre := model$coef[1] + model$coef[2] * year]
    }
}

for(i in 1:length(constituencies)) {
    message(i)
    make_pred_cons(constituencies[i])
    make_pred_cons_pre(constituencies[i])

}

# constructing leads of shock variable
df_s <- df_s[order(constituency_id, year)]
df_s[!(year %in% c(1880, 1890)), `:=`(shock_diff00w_lead2 = shift(shock_diff00w, -2),
    shock_diff00w_lead3 = shift(shock_diff00w, -3)),
    by = .(constituency_id)]



r_trend1s <- felm(cons_share ~ shock_diff00w + cons_share_pred |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900 & unopposed == 0])

r_trend2s <- felm(cons_share ~ shock_diff00w + cons_share_pred +
    as.factor(year) : const_frac_secondary|
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900 & unopposed == 0])

r_trend3s <- felm(cons_share ~ shock_diff00w + cons_share_pred_pre |
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900 & unopposed == 0])

r_trend4s <- felm(cons_share ~ shock_diff00w + cons_share_pred_pre +
    as.factor(year) : const_frac_secondary|
    constituency_id + year | 0 | county_name,
    df_s[year >= 1900 & unopposed == 0])

r_trend5s <- felm(cons_share ~ shock_diff00w_lead2  |
    constituency_id + year | 0 | county_name,
    df_s[year <= 1895 & unopposed == 0])

r_trend6s <- felm(cons_share ~ shock_diff00w_lead2 +
    as.factor(year) : const_frac_secondary |
    constituency_id + year | 0 | county_name,
    df_s[year <= 1895 & unopposed == 0])

r_trend7s <- felm(cons_share ~ shock_diff00w_lead3  |
    constituency_id + year | 0 | county_name,
    df_s[year <= 1895 & unopposed == 0])

r_trend8s <- felm(cons_share ~ shock_diff00w_lead3 +
    as.factor(year) : const_frac_secondary |
    constituency_id + year | 0 | county_name,
    df_s[year <= 1895 & unopposed == 0])


models <- list(r_trend1s, r_trend2s, r_trend3s, r_trend4s, r_trend5s, r_trend6s,
    r_trend7s, r_trend8s)
table_votes_trend_s <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$",
        "Constituency Conservative \\\\trend (1885--1910)",
        "Constituency Conservative \\\\trend (1885--1900)",
        "$\\Delta\\textrm{IPW}_{t+2}$",
        "$\\Delta\\textrm{IPW}_{t+3}$"),
    dep.var.labels.include = F,
    #keep = "shock",
    omit = "factor",
    dep.var.caption = "",
    add.lines = list(c("Years", "1900--1910", "1900--1910", "1900--1910",
        "1900--1910", "1885--1895", "1885--1895", "1885--1895", "1885--1895"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.7\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party. Models (1)--(4) re-estimate the
        main voting result, for the 1900--1910 period, adding controls for
        Conservative vote share as predicted by constituency-specific time
        trends, based on the 1885--1910 period for models (1) and (2), and based
        on the 1885--1900 period for (3) and (4). Models (5)--(8) test for
        differential trends in Conservative voting prior to the acceleration of
        German imports after 1895, regressing Conservative vote share 1885--1895
        on import penetration 1892--1906 and 1895--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("Checks for pre-trends"),
    label = "table_votes_trend_cl",
    notes.align = "l")
cat(table_votes_trend_s, sep = '\n',
    file = "Tables/table_votes_trend_cl.tex")

table_votes_trend_s_full <- stargazer(models,
    order = c("shock_diff00w", "cons_share_pred", "cons_share_pred_pre",
        "shock_diff00w_lead2", "shock_diff00w_lead3"),
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{IPW}_{t+2}$",
        "$\\Delta\\textrm{IPW}_{t+3}$",
        "Constituency Conservative \\\\trend (1885--1910)",
        "Constituency Conservative \\\\trend (1885--1900)"),
    dep.var.labels.include = F,
    #keep = "shock",
    #omit = "factor",
    dep.var.caption = "",
    add.lines = list(c("Years", "1900--1910", "1900--1910", "1900--1910",
        "1900--1910", "1885--1895", "1885--1895", "1885--1895", "1885--1895"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.7\\textwidth}{
        Constituency-level fixed effects regression, dependent variable is share
        of the vote for the Conservative Party. Models (1)--(4) re-estimate the
        main voting result, for the 1900--1910 period, adding controls for
        Conservative vote share as predicted by constituency-specific time
        trends, based on the 1885--1910 period for models (1) and (2), and based
        on the 1885--1900 period for (3) and (4). Models (5)--(8) test for
        differential trends in Conservative voting prior to the acceleration of
        German imports after 1895, regressing Conservative vote share 1885--1895
        on import penetration 1892--1906 and 1895--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-14"),
    label = "table_votes_trend_cl",
    notes.align = "l")
cat(table_votes_trend_s_full, sep = '\n',
    file = "Tables with full output/table_a14.tex")



# A-37: pre-trend checks with exposure-robust SEs
df_trend <- df_s[, .(constituency_id, year, cons_share_pred, cons_share_pred_pre)]
df_trend <- merge(df, df_trend, by = c("constituency_id", "year"), all.x = T)
unique(df$year)
df_trend2 <- df[year %in% c(1885, 1886, 1892, 1895, 1900, 1906, 1910, 1911)]
setorder(df_trend2, year)
df_trend2[, `:=`(diff_percap00w_lead2 = shift(diff_percap00w, -2),
    shock_diff00w_lead2 = shift(shock_diff00w, -2),
    diff_percap00w_lead3 = shift(diff_percap00w, -3),
    shock_diff00w_lead3 = shift(shock_diff00w, -3)),
    by = .(constituency_id, category_uk)]


r_trend1 <- borusyak_reg(data = df_trend[unopposed == 0 & year >= 1900 & !is.na(cons_share_pred)],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("cons_share_pred"),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_trend2 <- borusyak_reg(data = df_trend[unopposed == 0 & year >= 1900 & !is.na(cons_share_pred)],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("cons_share_pred", sec_controls),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_trend3 <- borusyak_reg(data = df_trend[unopposed == 0 & year >= 1900 &
    !is.na(cons_share_pred_pre)],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("cons_share_pred_pre"),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_trend4 <- borusyak_reg(data = df_trend[unopposed == 0 & year >= 1900 &
        !is.na(cons_share_pred_pre)],
    y_var = "cons_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    control_vars = c("cons_share_pred_pre", sec_controls),
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_trend5 <- borusyak_reg(data = df_trend2[unopposed == 0 & between(year, 1885, 1895)],
    y_var = "cons_share",
    x_var = "shock_diff00w_lead2",
    shock_var = "diff_percap00w_lead2",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
r_trend6 <- borusyak_reg(data = df_trend2[unopposed == 0 & between(year, 1885, 1895)],
    y_var = "cons_share",
    x_var = "shock_diff00w_lead2",
    shock_var = "diff_percap00w_lead2",
    share_var = "const_frac_category",
    control_vars = sec_controls,
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

r_trend7 <- borusyak_reg(data = df_trend2[unopposed == 0 & between(year, 1885, 1895)],
    y_var = "cons_share",
    x_var = "shock_diff00w_lead3",
    shock_var = "diff_percap00w_lead3",
    share_var = "const_frac_category",
    #control_vars = sec_controls,
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))



r_trend8 <- borusyak_reg(data = df_trend2[unopposed == 0 & between(year, 1885, 1895)],
    y_var = "cons_share",
    x_var = "shock_diff00w_lead3",
    shock_var = "diff_percap00w_lead3",
    share_var = "const_frac_category",
    control_vars = sec_controls,
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))

models <- list(r_trend1, r_trend2, r_trend3, r_trend4, r_trend5, r_trend6,
    r_trend7, r_trend8)
f_stats <- sapply(models, get_f_stat)

table_votes_trend <- stargazer(models,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}$"),
    dep.var.labels.include = F,
    keep = "x",
    #omit = "factor",
    dep.var.caption = "",
    add.lines = list(c("Years", "1900--1910", "1900--1910", "1900--1910",
        "1900--1910", "1885--1895", "1885--1895", "1885--1895", "1885--1895"),
        c("Shock", "$\\Delta\\textrm{IPW}_{1900}$", "$\\Delta\\textrm{IPW}_{1900}$",
            "$\\Delta\\textrm{IPW}_{1900}$", "$\\Delta\\textrm{IPW}_{1900}$",
            "$\\Delta\\textrm{IPW}_{t+2}$", "$\\Delta\\textrm{IPW}_{t+2}$",
            "$\\Delta\\textrm{IPW}_{t+3}$", "$\\Delta\\textrm{IPW}_{t+3}$"),
        c("Constituency \\\\time trend", "1885--1910", "1885--1910", "1885--1900",
            "1885--1900"),
        c("Initial MF x election", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-state", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.7\\textwidth}{
        This table replicates the results of Table \\ref{table_votes_trend_cl}
        using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022)
        Constituency-level variables aggregated up to the industry level,
        dependent variable is share
        of the vote for the Conservative Party. Models (1)--(4) re-estimate the
        main voting result, for the 1900--1910 period, adding controls for
        Conservative vote share as predicted by constituency-specific time
        trends, based on the 1885--1910 period for models (1) and (2), and based
        on the 1885--1900 period for (3) and (4). Models (5)--(8) test for
        differential trends in Conservative voting prior to the acceleration of
        German imports after 1895, regressing Conservative vote share 1885--1895
        on import penetration 1892--1906 and 1895--1910.
        All models include constituency and election fixed effects, even numbers
        add manufacturing employment in 1880 interacted with
        election dummies.
        Standard errors clustered by industry in parentheses.}",
    title = c("Checks for pre-trends, with exposure-robust standard errors"),
    label = "table_votes_trend_ex",
    notes.align = "l")
cat(table_votes_trend, sep = '\n',
    file = "Tables/table_votes_trend_ex.tex")
rm(df_trend, df_trend2)

# A-15: robustness checkf with Chaisemartin-d'Haultfoeuille estimator

# aggregating 1910 elections
df_s[, `:=`(const_frac_secondary_1885 = const_frac_secondary * (year == 1885),
    const_frac_secondary_1886 = const_frac_secondary * (year == 1886),
    const_frac_secondary_1892 = const_frac_secondary * (year == 1892),
    const_frac_secondary_1895 = const_frac_secondary * (year == 1895),
    const_frac_secondary_1906 = const_frac_secondary * (year == 1906),
    const_frac_secondary_1910 = const_frac_secondary * (year == 1910),
    const_frac_secondary_1911 = const_frac_secondary * (year == 1911))]

df_to_bs <- df_s[unopposed == 0, .(constituency_id, county_name, year,
    shock_diff_w, shock_diff00w, cons_share,
    const_frac_secondary_1885,
    const_frac_secondary_1886, const_frac_secondary_1892,
    const_frac_secondary_1895,
    const_frac_secondary_1906, const_frac_secondary_1910,
    const_frac_secondary_1911)]

df_to_bs[year == 1911, `:=`(year = 1910,
    const_frac_secondary_1910 = const_frac_secondary_1911)]

df_to_bs <- df_to_bs[, .(cons_share = mean(cons_share, na.rm = TRUE)),
    by = .(constituency_id, county_name, year, shock_diff_w,
        shock_diff00w, const_frac_secondary_1885,
        const_frac_secondary_1886, const_frac_secondary_1892,
        const_frac_secondary_1895,
        const_frac_secondary_1906, const_frac_secondary_1910)]

# rounding shock to the nearest 0.5
df_to_bs[, shock_round := round(shock_diff00w * 2, 0) / 2]
dim(df_to_bs)


r_ch1 <- felm(cons_share ~ shock_round
    | constituency_id + year | 0 | county_name,
    df_to_bs[year >= 1900])

r_ch2 <- felm(cons_share ~ shock_round + const_frac_secondary_1906 +
    const_frac_secondary_1910
    | constituency_id + year | 0 | county_name,
    df_to_bs[year >= 1900])




r_ch3 <- did_multiplegt(df = df_to_bs[year >= 1900 &
    !is.na(shock_round) &
    !is.na(cons_share) & !is.na(county_name)],
    Y = "cons_share",
    G = "constituency_id",
    T = "year",
    D = "shock_round"
)

str(r_ch3)

county_weights <- df_to_bs[year >= 1900 & !is.na(shock_diff00w), .(n = .N), by = .(county_name)]
county_weights <- county_weights[order(county_name)]
df_to_bs <- df_to_bs[order(constituency_id, year)]


# function to bootstrap
# Note: bootstrapping the CH estimator generally takes around an hour

RNGkind("L'Ecuyer-CMRG")
bootstrap_r_ch3 <- function() {

    county_sample <- county_weights$county_name[
        sample(x = nrow(county_weights),
            size = nrow(county_weights),
            replace = TRUE)]

    df_sample <- lapply(county_sample, function(x) {
        return(df_to_bs[county_name == x & year >= 1900
            ])
    }) %>% rbindlist()
    did_bs <- did_multiplegt(df = df_sample,
        Y = "cons_share",
        G = "constituency_id",
        T = "year",
        D = "shock_round"
    )
    return(did_bs$effect)
}
set.seed(123)

bs_r_ch3 <- mclapply(1:1000, function(x) {
    message(x)
    bs_effect <- bootstrap_r_ch3()
    return(bs_effect)
}, mc.cores = 3, mc.set.seed = TRUE)
bs_r_ch3 <- simplify2array(bs_r_ch3)

confint_r_ch3 <- quantile(bs_r_ch3, c(0.025, 0.975))
print(confint_r_ch3)

r_ch4 <- did_multiplegt(df = df_to_bs[year >= 1900 &
    !is.na(shock_round) &
    !is.na(cons_share) & !is.na(county_name)],
    Y = "cons_share",
    G = "constituency_id",
    T = "year",
    D = "shock_round",
    controls = c("const_frac_secondary_1906", "const_frac_secondary_1910")
)


bootstrap_r_ch4 <- function() {

    county_sample <- county_weights$county_name[
        sample(x = nrow(county_weights),
            size = nrow(county_weights),
            #prob = county_weights$prob,
            replace = TRUE)]

    df_sample <- lapply(county_sample, function(x) {
        return(df_to_bs[county_name == x & year >= 1900
            ])
    }) %>% rbindlist()
    did_bs <- did_multiplegt(df = df_sample,
        Y = "cons_share",
        G = "constituency_id",
        T = "year",
        D = "shock_round",
        controls = c(
            "const_frac_secondary_1906",
            "const_frac_secondary_1910"
        )
    )
    return(did_bs$effect)
}
set.seed(123)
bs_r_ch4 <- mclapply(1:1000, function(x) {
    message(x)
    bs_effect <- bootstrap_r_ch4()
    return(bs_effect)
}, mc.cores = 3, mc.set.seed = TRUE)
bs_r_ch4 <- simplify2array(bs_r_ch4)

confint_r_ch3 <- quantile(bs_r_ch3, c(0.025, 0.975))
confint_r_ch4 <- quantile(bs_r_ch4, c(0.025, 0.975))

confint_r_ch4

mod_ch1 <- createTexreg(coef.names = rownames(r_ch1$coefficients),
    coef = r_ch1$coefficients[, 1],
    se = r_ch1$cse,
    pvalues = r_ch1$cpval,
    gof.names = c("N"),
    gof = c(r_ch1$N),
    gof.decimal = c(FALSE)
    )
mod_ch2 <- createTexreg(coef.names = rownames(r_ch2$coefficients),
    coef = r_ch2$coefficients[, 1],
    se = r_ch2$cse,
    pvalues = r_ch2$cpval,
    gof.names = c("N"),
    gof = c(r_ch2$N),
    gof.decimal = c(FALSE)
    )
mod_ch3 <- createTexreg(coef.names = "shock_round",
    coef = r_ch3$effect,
    ci.low = confint_r_ch3[1],
    ci.up = confint_r_ch3[2],
    pvalues = mean(bs_r_ch3 > 0) / 2,
    gof.names = c("N", "N switchers"),
    gof = c(r_ch3$N_effect, r_ch3$N_switchers_effect),
    gof.decimal = c(FALSE, FALSE)
)
mod_ch4 <- createTexreg(coef.names = "shock_round",
    coef = r_ch4$effect,
    ci.low = confint_r_ch4[1],
    ci.up = confint_r_ch4[2],
    pvalues = mean(bs_r_ch4 > 0) / 2,
    gof.names = c("N", "N switchers"),
    gof = c(r_ch4$N_effect, r_ch4$N_switchers_effect),
    gof.decimal = c(FALSE, FALSE)
)



texreg(list(mod_ch1, mod_ch2, mod_ch3, mod_ch4),
    file = "Tables/table_ch_robust.tex",
    label = "table_ch_robust",
    omit.coef = "const",
    custom.coef.names = "$\\Delta\\textrm{IPW}_{1900}$ (rounded)",
    custom.model.names = c("(1)", "(2)", "(3)", "(4)"),
    caption = "Robustness of post-1900 voting results to
        Chaisemartin-D'Haultfoeuille estimator",
    custom.gof.rows = list("Initial Mf x year" = c("", "x", "", "x")),
    custom.header = list("TWFE estimator" = 1:2,
        "CH estimator" = 3:4),
    custom.note = paste0("\\item This table shows the results of regressions of ",
        "Conservative vote share, 1900--1910 on the change in imports per worker. ",
        "Models (1) and (2) use the conventional ",
        "two-way fixed effects estimator used throughout the article. Models (3) ",
        "and (4) use the estimator proposed by Chaisemartin and D'Haultfoeuille, ",
        "which corrects for negative weights. This estimator directly compares units ",
        "which changed treatment status from one period to the next against units ",
        "which did not. In order to use this estimator, we round our ",
        "$\\Delta\\textrm{IPW}$ measure to the nearest 0.5, and average the ",
        "dependent variable over the two 1910 elections (for which the treatment ",
        "is unchanged). All models control for constituency and year fixed effects, ",
        "and (2) and (4) control for initial manufacturing interacted with year ",
        "fixed effects. For models (1) and (2), standard errors clustered by ",
        "county are shown in parentheses, for (3) and (4) we cluster bootstrap ",
        "at the county level and report 95\\% confidence intervals. ",
        " %stars"),    #ci.force = TRUE,
    stars = c(0.05),
    #star.symbol = "\\*\\*",
    #ci.test = NA,
    digits = 3,
    booktabs = T,
    threeparttable = T,
    #sideways = T,
    use.packages = F
    )

# Table A-16 Incumbency effects
ri_1_1s <- felm(incumbent_mp_share ~ shock_diff_w | constituency_id + year |
    0 | county_name,
    data = df_s[unopposed == 0])
ri_1_2s <- felm(incumbent_mp_share ~ shock_diff00w | constituency_id + year |
    0 | county_name,
    data = df_s[year >= 1900 & unopposed == 0])
ri_2_1s <- felm(incumbent_party_share ~ shock_diff_w | constituency_id + year |
    0 | county_name,
    data = df_s[unopposed == 0])
ri_2_2s <- felm(incumbent_party_share ~ shock_diff00w | constituency_id + year |
    0 | county_name,
    data = df_s[year >= 1900 & unopposed == 0])
ri_3_1s <- felm(fd_incum ~ fd_shock_w |  year |
    0 | county_name,
    data = df_s[!is.na(fd_incum) & year >= 1886])
ri_3_2s <- felm(fd_incum ~ fd_shock_w |  year |
    0 | county_name,
    data = df_s[!is.na(fd_incum) & year >1900])
table_incum_s <- stargazer(ri_1_1s, ri_1_2s, ri_2_1s, ri_2_2s, ri_3_1s, ri_3_2s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$",
        "$\\Delta\\textrm{IPW}_{t}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = "Constant",
    column.labels = c("MP", "Local Party", "National Party"),
    column.separate = c(2, 2, 2),
    add.lines = list(c("Years", "All", "1900-1910", "All",  "1900-1910",
        "All",  "1900-1910")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Constituency-level regressions,
        (1)--(4) are estimated in levels and include constituency and year fixed effects,
        (5) and (6) in stacked first-differences, and include year fixed effects.
        For (1) and (2) the dependent variable is the share of the vote
        won by incumbent MPs, for (3) and (4), the share of the vote won by incumbent
        parties at the local level, for (5) and (6), the change in voteshare by the
        nationally-incumbent party. Standard errors clustered by county in
        parentheses.}",
    title = c("Effects of import competition on incumbency"),
    label = "table_incumbency_cl",
    notes.align = "l")
cat(table_incum_s, sep = '\n', file = "Tables/table_incumbency_cl.tex")
rm(ri_1_1s, ri_1_2s, ri_2_1s, ri_2_2s, ri_3_1s, ri_3_2s)

# A-38: incumbency with exposure-robust SEs
ri_1_1 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "incumbent_mp_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
ri_1_2 <- borusyak_reg(data = df[year >= 1900 & unopposed == 0],
    y_var = "incumbent_mp_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
ri_2_1 <- borusyak_reg(data = df[unopposed == 0],
    y_var = "incumbent_party_share",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
ri_2_2 <- borusyak_reg(data = df[year >= 1900 & unopposed == 0],
    y_var = "incumbent_party_share",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("constituency_id", "year"))
ri_3_1 <- borusyak_reg(data = df[year >= 1886 & !is.na(fd_incum)],
    y_var = "fd_incum",
    x_var = "fd_shock_w",
    shock_var = "fd_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"))
ri_3_2 <- borusyak_reg(data = df[year > 1900 & !is.na(fd_incum)],
    y_var = "fd_incum",
    x_var = "fd_shock_w",
    shock_var = "fd_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("year"))
table_incum <- stargazer(ri_1_1, ri_1_2, ri_2_1, ri_2_2, ri_3_1, ri_3_2,
    covariate.labels = c("$\\Delta$IPW"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = "Constant",
    column.labels = c("MP", "Local Party", "National Party"),
    column.separate = c(2, 2, 2),
    add.lines = list(c("Years", "All", "1900-1910", "All",  "1900-1910",
        "All",  "1900-1910"),
        c("First stage F-stat", get_f_stat(ri_1_1), get_f_stat(ri_1_2),
        get_f_stat(ri_2_1), get_f_stat(ri_2_2), get_f_stat(ri_3_1), get_f_stat(ri_3_2))),
    keep.stat = c("n"),
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_incumbency_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Constituency-level variables aggregated to the industry level,
        (1)--(4) are estimated in levels and include constituency and year fixed effects,
        (5) and (6) in stacked first-differences, and include year fixed effects.
        For (1) and (2) the dependent variable is the share of the vote
        won by incumbent MPs, for (3) and (4), the share of the vote won by incumbent
        parties at the local level, for (5) and (6), the change in voteshare by the
        nationally-incumbent party. Standard errors clustered by industry in
        parentheses.}",
    title = c("Effects of import competition on incumbency,
        exposure-robust standard errors"),
    label = "table_incumbency_ex",
    notes.align = "l")
cat(table_incum, sep = '\n', file = "Tables/table_incumbency_ex.tex")
rm(ri_1_1, ri_1_2, ri_2_1, ri_2_2, ri_3_1, ri_3_2)

#%% Newspaper regressions
# Table 4: effects of IPW on references to trade
rt_ts1s <- shock_term_s(df_news[term %in% "import"])
rt_ts2s <- shock_term00s(df_news[term %in% "import" & year >= 1900])
rt_ts3s <- shock_term_s(df_news[term %in% "trade"])
rt_ts4s <- shock_term00s(df_news[term %in% "trade" & year >= 1900])



models <- list(rt_ts1s[[1]], rt_ts1s[[2]], rt_ts2s[[1]], rt_ts2s[[2]],
    rt_ts3s[[1]], rt_ts3s[[2]], rt_ts4s[[1]], rt_ts4s[[2]])
table_news_trade_short_s <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``import''",  "``trade''"),
    column.separate = c(4, 4),
    add.lines = list(c("Years", "All", "All", "1900--1910",
        "1900--1910", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on newspaper references to trade"),
    label = "table_news_trade_short_cl",
    notes.align = "l")
cat(table_news_trade_short_s, sep = '\n',
    file = "Tables/table_news_trade_short_cl.tex")

table_news_trade_short_s_full <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    order = c("shock_diff_w", "shock_diff00w"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``import''",  "``trade''"),
    column.separate = c(4, 4),
    add.lines = list(c("Years", "All", "All", "1900--1910",
        "1900--1910", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table 4"),
    label = "table_news_trade_short_cl",
    notes.align = "l")
cat(table_news_trade_short_s_full, sep = '\n',
    file = "Tables with full output/table_4.tex")


# Table A-24
rt_ts1 <- shock_term(df_news[term %in% "import"])
rt_ts2 <- shock_term00(df_news[term %in% "import" & year >= 1900])
rt_ts3 <- shock_term(df_news[term %in% "trade"])
rt_ts4 <- shock_term00(df_news[term %in% "trade" & year >= 1900])
models <- list(rt_ts1[[1]], rt_ts1[[2]], rt_ts2[[1]], rt_ts2[[2]],
    rt_ts3[[1]], rt_ts3[[2]], rt_ts4[[1]], rt_ts4[[2]])
f_stats <- sapply(models, get_f_stat)
table_news_trade_short <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``import''",  "``trade''"),
    column.separate = c(4, 4),
    add.lines = list(c("Years", "All", "All", "1900--1910",
        "1900--1910", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{This table replicates the results of
        Table \\ref{table_news_trade_short_cl} using the aggregation and
        standard error calculation methods recommended by Borusyak, Hull, and
        Jaravel (2022).
        Newspaper-level variables aggregated up to the industry level.
        Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on newspaper references to trade,
        exposure-robust standard errors"),
    label = "table_news_trade_short_ex",
    notes.align = "l")
cat(table_news_trade_short, sep = '\n',
    file = "Tables/table_news_trade_short_ex.tex")

rm(rt_1, rt_2, rt_3, rt_4, rt_5, rt_1s, rt_2s, rt_3s, rt_4s, rt_5s,
    rt_t1s, rt_t2s, rt_t3s, rt_t4s, rt_t1, rt_t2, rt_t3, rt_t4,
    rt_1r, rt_2r, rt_3r, rt_4r, rt_5r, rt_ts1, rt_ts2, rt_ts3, rt_ts4,
    rt_ts1s, rt_ts2s, rt_ts3s, rt_ts4s)


# Table 6: relative usage of terms related to unemployment and vagrancy
unique(df_news[year == 1880, term])
# creating standardized measure of references to unemployment vs vagrancy
df_unemp_terms <- df_news[term %in% c("vagrant", "vagrancy", "pauper",
    "pauperism", "unemployed", "unemployment", "employment"),
    .(unemp = sum(per_issue * term %in% c("unemployed", "unemployment", "employment")),
    paup_vag = sum(per_issue * term %in% c("vagrant", "pauper", "vagrancy", "pauperism"))),
    by = .(title, constituency_group, year)]

df_unemp_terms[, unemp_diff := unemp - paup_vag]
df_unemp_terms[, unemp_diff_sd := (unemp_diff - mean(unemp_diff)) / sd(unemp_diff)]

df_unemp_terms_s <- merge(df_unemp_terms, df_group_s, by = c("constituency_group", "year"))
df_unemp_terms <- merge(df_unemp_terms, df_group, by = c("constituency_group", "year"))
colnames(df_unemp_terms_s)



r_unemp1s <- felm(unemp_diff_sd ~ shock_diff_w | title + year | 0 | county_name,
    df_unemp_terms_s[])

r_unemp2s <- felm(unemp_diff_sd ~ shock_diff_w + const_frac_secondary : as.factor(year)
    | title + year | 0 | county_name,
    df_unemp_terms_s[])

r_unemp3s <- felm(unemp_diff_sd ~ shock_diff00w | title + year | 0 | county_name,
    df_unemp_terms_s[year >= 1900])

r_unemp4s <- felm(unemp_diff_sd ~ shock_diff00w + const_frac_secondary :
    as.factor(year)| title + year | 0 | county_name,
    df_unemp_terms_s[year >= 1900])

table_unemp_terms_s <- stargazer(r_unemp1s, r_unemp2s, r_unemp3s, r_unemp4s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    add.lines = list(c("Years", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        Newspaper-level regressions. Dependent variable is the number of references to
        ``unemployed,'' ``unemployment,'' and ``employment,'' minus the number of
        references to ``vagrants,'' ``vagrancy,'' ``pauper,'' and ``pauperism,''
        standardized. All models include newspaper and year fixed effects.
        For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on newspaper references to
        unemployment, vagrancy, and pauperism"),
    label = "table_unemp_terms_cl",
    notes.align = "l")
cat(table_unemp_terms_s, sep = '\n', file = "Tables/table_unemp_terms_cl.tex")

table_unemp_terms_s_full <- stargazer(r_unemp1s, r_unemp2s, r_unemp3s, r_unemp4s,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$",
        "$\\Delta\\textrm{IPW}_{1900}$"),
    order = c("shock_diff_w", "shock_diff00w"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    add.lines = list(c("Years", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.3\\textwidth}{
        Newspaper-level regressions. Dependent variable is the number of references to
        ``unemployed,'' ``unemployment,'' and ``employment,'' minus the number of
        references to ``vagrants,'' ``vagrancy,'' ``pauper,'' and ``pauperism,''
        standardized. All models include newspaper and year fixed effects.
        For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table 6"),
    label = "table_unemp_terms_cl",
    notes.align = "l")
cat(table_unemp_terms_s_full, sep = '\n', file = "Tables with full output/table_6.tex")

rm(r_unemp1s, r_unemp2s, r_unemp3s, r_unemp4s)


# Table A-26: replicating this with exposure-robust SEs

r_unemp1 <- borusyak_reg(data = df_unemp_terms[],
    y_var = "unemp_diff_sd",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("title", "year"))
r_unemp2 <- borusyak_reg(data = df_unemp_terms[],
    y_var = "unemp_diff_sd",
    x_var = "shock_diff_w",
    shock_var = "diff_percap_w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("title", "year"),
    control_vars = sec_controls_s)
r_unemp3 <- borusyak_reg(data = df_unemp_terms[year >= 1900],
    y_var = "unemp_diff_sd",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("title", "year"))
r_unemp4 <- borusyak_reg(data = df_unemp_terms[year >= 1900],
    y_var = "unemp_diff_sd",
    x_var = "shock_diff00w",
    shock_var = "diff_percap00w",
    share_var = "const_frac_category",
    group_vars = c("category_uk", "year"), # entity, year
    fe_vars = c("title", "year"),
    control_vars = sec_controls_s)



table_unemp_terms <- stargazer(r_unemp1, r_unemp2, r_unemp3, r_unemp4,
    covariate.labels = c("$\\Delta\\textrm{IPW}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    add.lines = list(c("Years", "All", "All", "1900--1910", "1900--1910"),
        c("Initial Mf x year", "", "x", "", "x"),
        c("First stage F-state", get_f_stat(r_unemp1), get_f_stat(r_unemp2),
            get_f_stat(r_unemp3), get_f_stat(r_unemp4))),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        This table replicates the results of Table \\ref{table_unemp_terms_cl}
        using the aggregation and standard error calculation methods
        recommended by Borusyak, Hull, and Jaravel (2022). Newspaper-level
        variables aggregated to the industry level.
        Dependent variable is the number of references to
        ``unemployed,'' ``unemployment,'' and ``employment,'' minus the number of
        references to ``vagrants,'' ``vagrancy,'' ``pauper,'' and ``pauperism,''
        standardized. All models include newspaper and year fixed effects.
        For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on newspaper references to
        unemployment, vagrancy, and pauperism, exposure-robust standard errors"),
    label = "table_unemp_terms_ex",
    notes.align = "l")
cat(table_unemp_terms, sep = '\n', file = "Tables/table_unemp_terms_ex.tex")
rm(r_unemp1, r_unemp2, r_unemp3, r_unemp4)

# Table A-17: effect of IPW on terms used in Beveridge's analysis of unemployment
beveridge_terms <- c("unemployed", "unemployment", "industrial", "exchange",
    "table", "fluctuation", "demand", "depression", "trades", "reserve",
    "percentage", "organisation", "situation", "cyclical", "skilled", "dock",
    "note", "seasonal", "unskilled", "production")
rb_1s <- shock_term_s(df_news[term %in% beveridge_terms])
rb_2s <- shock_term_s(df_news[term %in% "fluctuation"])
rb_3s <- shock_term_s(df_news[term %in% "depression"])

rb_4s <- shock_term_s(df_news[term %in% c("unemployment", "employment")])
rb_5s <- shock_term_s(df_news[term %in% c("exchange")])
table_news_beveridge_s <- stargazer(rb_1s[[1]], rb_1s[[2]], rb_2s[[1]], rb_2s[[2]],
    rb_3s[[1]], rb_3s[[2]], rb_4s[[1]], rb_4s[[2]], rb_5s[[1]], rb_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("Beveridge terms", "``fluctuation''", "``depression''",
        "``(un)employment''", "``exchange''"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``Beveridge terms''
        refers to terms overused in Beveridge's \\emph{Unemployment: A Problem
        of Industry}, relative to other contemporary writings supportive of the
        existing Poor Law system. Terms were selected using the $\\chi^2$ test
        statistic proposed by Gentzkow and Shapiro (2010). The terms in question
        are ``unemployed,'' ``unemployment,'' ``industrial,'' ``exchange,''
        ``table,'' ``fluctuation,'' ``demand,'' ``depression,'' ``trades,''
        ``reserve,'' ``percentage,'' ``organisation,'' ``situation,''
        ``cyclical,'' ``skilled,'' ``dock,'' ``note,'' ``seasonal,''
        ``unskilled,'' and ``production.'' Standard errors clustered by county in
        parentheses.}",
    title = c("Effects of import competition on newspaper references terms
        overused in Beveridge's analysis of unemployment"),
    label = "table_news_beveridge_cl",
    notes.align = "l")
cat(table_news_beveridge_s, sep = '\n', file = "Tables/table_news_beveridge_cl.tex")

table_news_beveridge_s_full <- stargazer(rb_1s[[1]], rb_1s[[2]], rb_2s[[1]], rb_2s[[2]],
    rb_3s[[1]], rb_3s[[2]], rb_4s[[1]], rb_4s[[2]], rb_5s[[1]], rb_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("Beveridge terms", "``fluctuation''", "``depression''",
        "``(un)employment''", "``exchange''"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.7\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``Beveridge terms''
        refers to terms overused in Beveridge's \\emph{Unemployment: A Problem
        of Industry}, relative to other contemporary writings supportive of the
        existing Poor Law system. Terms were selected using the $\\chi^2$ test
        statistic proposed by Gentzkow and Shapiro (2010). The terms in question
        are ``unemployed,'' ``unemployment,'' ``industrial,'' ``exchange,''
        ``table,'' ``fluctuation,'' ``demand,'' ``depression,'' ``trades,''
        ``reserve,'' ``percentage,'' ``organisation,'' ``situation,''
        ``cyclical,'' ``skilled,'' ``dock,'' ``note,'' ``seasonal,''
        ``unskilled,'' and ``production.'' Standard errors clustered by county in
        parentheses.}",
    title = c("Full regression output for Table A-17"),
    label = "table_news_beveridge_cl",
    notes.align = "l")
cat(table_news_beveridge_s_full, sep = '\n',
    file = "Tables with full output/Table_a17.tex")



# A-39: Beveridge terms with exposure-robust SEs
rm(rt_1s, rt_2s, rt_3s, rt_4s, rt_5s)

rb_1 <- shock_term(df_news[term %in% beveridge_terms])
rb_2 <- shock_term(df_news[term %in% "fluctuation"])
rb_3 <- shock_term(df_news[term %in% "depression"])
rb_4 <- shock_term(df_news[term %in% c("unemployment", "employment")])
rb_5 <- shock_term(df_news[term %in% c("exchange")])
models <- list(rb_1[[1]], rb_1[[2]], rb_2[[1]], rb_2[[2]],
    rb_3[[1]], rb_3[[2]], rb_4[[1]], rb_4[[2]], rb_5[[1]], rb_5[[2]])
f_stats <- sapply(models, get_f_stat)
f_stats

table_news_beveridge <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("Beveridge terms", "``fluctuation''", "``depression''",
        "``(un)employment''", "``exchange''"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.8\\textwidth}{
        This table replicates the results of Table
        \\ref{table_news_beveridge_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Newspaper-level variables aggregated up to the
        industry level. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``Beveridge terms''
        refers to terms overused in Beveridge's \\emph{Unemployment: A Problem
        of Industry}, relative to other contemporary writings supportive of the
        existing Poor Law system. Terms were selected using the $\\chi^2$ test
        statistic proposed by Gentzkow and Shapiro (2010). The terms in question
        are ``unemployed,'' ``unemployment,'' ``industrial,'' ``exchange,''
        ``table,'' ``fluctuation,'' ``demand,'' ``depression,'' ``trades,''
        ``reserve,'' ``percentage,'' ``organisation,'' ``situation,''
        ``cyclical,'' ``skilled,'' ``dock,'' ``note,'' ``seasonal,''
        ``unskilled,'' and ``production.''
        Standard errors clustered by industry in
        parentheses.}",
    title = c("Effects of import competition on newspaper references terms
        overused in Beveridge's analysis of unemployment,
        exposure-robust standard errors"),
    label = "table_news_beveridge_ex",
    notes.align = "l")
cat(table_news_beveridge, sep = '\n',
    file = "Tables/table_news_beveridge_ex.tex")
rm(rb_1, rb_2, rb_3, rb_4, rb_5)

# Table A-18 effects of IPW on references to Germnay in newsappers
german_terms <- c("germany", "kaiser", "teuton", "prussia", "fatherland")
navy_terms <- c("navy", "naval", "dreadnought", "battleship", "fleet")
militarist_orgs <- c('"navy league"', '"national service league"')

rg1s <- shock_term_s(df_news[term %in% c("germany")])
rg2s <- shock_term_s(df_news[term %in% german_terms])
rg3s <- shock_term_s(df_news[term %in% navy_terms])
rg4s <- shock_term_s(df_news[term %in% militarist_orgs])
table_news_german_s <- stargazer(rg1s[[1]], rg1s[[2]], rg2s[[1]], rg2s[[2]],
    rg3s[[1]], rg3s[[2]], rg4s[[1]], rg4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``German terms'' are
        ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,'' and ``fatherland,''
        ``Navy terms'' are ``navy,'' ``naval,'' ``dreadnought,'' ``battleship,''
        and ``fleet,'' ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on newspaper references to Germany
        and the naval race"),
    label = "table_news_german_cl",
    notes.align = "l")
cat(table_news_german_s, sep = '\n', file = "Tables/table_news_german_cl.tex")

table_news_german_s_full <- stargazer(rg1s[[1]], rg1s[[2]], rg2s[[1]], rg2s[[2]],
    rg3s[[1]], rg3s[[2]], rg4s[[1]], rg4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.5\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``German terms'' are
        ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,'' and ``fatherland,''
        ``Navy terms'' are ``navy,'' ``naval,'' ``dreadnought,'' ``battleship,''
        and ``fleet,'' ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-18"),
    label = "table_news_german_cl",
    notes.align = "l")
cat(table_news_german_s_full, sep = '\n', file = "Tables with full output/table_a18.tex")


# A-40: effect of IPW on references to Germany, exposure-robust SEs
rg1 <- shock_term(df_news[term %in% c("germany")])
rg2 <- shock_term(df_news[term %in% german_terms])
rg3 <- shock_term(df_news[term %in% navy_terms])
rg4 <- shock_term(df_news[term %in% militarist_orgs])
models <- list(rg1[[1]], rg1[[2]], rg2[[1]], rg2[[2]],
    rg3[[1]], rg3[[2]], rg4[[1]], rg4[[2]])
f_stats <- sapply(models, get_f_stat)
table_news_german <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_news_german_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Newspaper-level data aggregated to the industry level.
        Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level. ``German terms'' are
        ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,'' and ``fatherland,''
        ``Navy terms'' are ``navy,'' ``naval,'' ``dreadnought,'' ``battleship,''
        and ``fleet,'' ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on newspaper references to Germany
        and the naval race, exposure-robust standard errors"),
    label = "table_news_german_ex",
    notes.align = "l")
cat(table_news_german, sep = '\n', file = "Tables/table_news_german_ex.tex")
rm(models, rg1, rg2, rg3, rg4, rg1s, rg2s, rg3s, rg4s)

# Table A-21: local news coverage of immigration
ran_1s <- shock_term_alien_s(df_news[term %in% c("immigrant")])
ran_2s <- shock_term_alien_s(df_news[term %in% c("alien")])
ran_3s <- shock_term_alien_s(df_news[term %in% c("jew")])
ran_4s <- shock_term_alien_s(df_news[term %in% c("foreigner")])
ran_5s <- shock_term_alien_s(df_news[
    term %in% c("immigrant", "alien", "jew", "foreigner")])
table_news_aliens_s <- stargazer(ran_1s[[1]], ran_1s[[2]], ran_2s[[1]],
    ran_2s[[2]], ran_3s[[1]], ran_3s[[2]], ran_4s[[1]], ran_4s[[2]],
    ran_5s[[1]], ran_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary", "const_frac_n_irish"),
    column.labels = c("``immigrant''", "``alien''", "``jew''", "``foreigner''", "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x", "x",
        "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        (9) and (10) use mentions
        of all four terms. Standard errors clustered by county in
        parentheses.}",
    title = c(
        "Effect of local trade shocks on newspaper coverage of immigration"),
    label = "table_news_aliens_cl",
    notes.align = "l"
)
cat(table_news_aliens_s, sep = '\n',
    file = "Tables/table_news_aliens_cl.tex")

table_news_aliens_s_full <- stargazer(ran_1s[[1]], ran_1s[[2]], ran_2s[[1]],
    ran_2s[[2]], ran_3s[[1]], ran_3s[[2]], ran_4s[[1]], ran_4s[[2]],
    ran_5s[[1]], ran_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary", "const_frac_n_irish"),
    column.labels = c("``immigrant''", "``alien''", "``jew''", "``foreigner''", "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x", "x",
        "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Newspaper-level regressions. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        (9) and (10) use mentions
        of all four terms. Standard errors clustered by county in
        parentheses.}",
    title = c(
        "Full regression output for Table A-21"),
    label = "table_news_aliens_cl",
    notes.align = "l"
)
cat(table_news_aliens_s_full, sep = '\n',
    file = "Tables with full output/table_a21.tex")
rm(ran_1s, ran_2s, ran_3s, ran_4s, ran_5s)

# Table A-43: local news coverage of immigration, exposure-robust SEs
ran_1 <- shock_term_alien(df_news[term %in% c("immigrant")])
ran_2 <- shock_term_alien(df_news[term %in% c("alien")])
ran_3 <- shock_term_alien(df_news[term %in% c("jew")])
ran_4 <- shock_term_alien(df_news[term %in% c("foreigner")])
ran_5 <- shock_term_alien(df_news[
    term %in% c("immigrant", "alien", "jew", "foreigner")])
models <- list(ran_1[[1]], ran_1[[2]], ran_2[[1]], ran_2[[2]],
    ran_3[[1]], ran_3[[2]], ran_4[[1]], ran_4[[2]], ran_5[[1]], ran_5[[2]])
f_stats <- sapply(models, get_f_stat)

table_news_aliens <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``immigrant'", "``alien''", "``jew''", "``foreigner''",
        "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x",
        "x", "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_news_aliens_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Newspaper-level variables aggregated to the
        industry level. Dependent variable is number of uses of
        specified term per newspaper issue, standardized. All models include
        newspaper and year fixed effects. For newspapers in cities, $\\Delta$IPW
        is calculated at the city-, not constituency-level.
        (9) and (10) use mentions
        of all four terms. Standard errors clustered by industry in
        parentheses.}",
    label = "table_news_aliens_ex",
    title = c("Effect of local trade shocks on newspaper coverage of immigration,
        exposure-robust standard errors"),
    notes.align = "l")
cat(table_news_aliens, sep = '\n', file = "Tables/table_news_aliens_ex.tex")
rm(models)
rm(ran_1, ran_2, ran_3, ran_4, ran_5)

#%% Regressions with manifesto data

# Table 5: effect of IPW on Liberal references to social reform
rl_1s <- manifesto_reg_lib_s(df_manifestos[term == "social reform"])
rl_2s <- manifesto_reg_lib_s(df_manifestos[term == "poor law"])
rl_3s <- manifesto_reg_lib_s(df_manifestos[term == "labour exchange"])
rl_4s <- manifesto_reg_lib_s(df_manifestos[term == "social reform|poor law|labour exchange"])

table_manif_lib_s <- stargazer(rl_1s[[1]], rl_1s[[2]], rl_2s[[1]], rl_2s[[2]],
    rl_3s[[1]], rl_3s[[2]], rl_4s[[1]], rl_4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``social reform''", "``poor law''", "``labour exchange''",
        "All"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Liberal
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by county in parentheses.}",
    title = c("Effect of local trade shocks on references to social reform in
        Liberal campaign manifestos"),
    label = "table_manif_lib_cl",
    notes.align = "l"
)
cat(table_manif_lib_s, sep = '\n', file = "Tables/table_manif_lib_cl.tex")
table_manif_lib_s_full <- stargazer(rl_1s[[1]], rl_1s[[2]], rl_2s[[1]], rl_2s[[2]],
    rl_3s[[1]], rl_3s[[2]], rl_4s[[1]], rl_4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``social reform''", "``poor law''", "``labour exchange''",
        "All"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.4\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Liberal
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table 5"),
    label = "table_manif_lib_cl",
    notes.align = "l"
)
cat(table_manif_lib_s_full, sep = '\n', file = "Tables with full output/table_5.tex")

rm(rl_1s, rl_2s, rl_3s, rl_4s)

# A-25: effects of IPW on Liberal social reform, with exposure-robust SEs
rl_1 <- manifesto_reg_lib(df_manifestos[term == "social reform"])
rl_2 <- manifesto_reg_lib(df_manifestos[term == "poor law"])
rl_3 <- manifesto_reg_lib(df_manifestos[term == "labour exchange"])
rl_4 <- manifesto_reg_lib(df_manifestos[term == "social reform|poor law|labour exchange"])
models <- list(rl_1[[1]], rl_1[[2]], rl_2[[1]], rl_2[[2]],
    rl_3[[1]], rl_3[[2]], rl_4[[1]], rl_4[[2]])
f_stats <- sapply(models, get_f_stat)

table_manif_lib <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``social reform''", "``poor law''", "``labour exchange''",
        "All"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_manif_lib_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022).
        Manifesto-level data aggregated to the industry level.
        Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Liberal
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effect of local trade shocks on references to social reform in
        Liberal campaign manifestos, exposure-robust standard errors"),
    label = "table_manif_lib_ex",
    notes.align = "l"
)
cat(table_manif_lib, sep = '\n', file = "Tables/table_manif_lib_ex.tex")

# Table A-19: effects of IPW on references to Germany and naval race in manifestos
rgm1s <- manifesto_reg_all_s(df_manifestos[term == "germany"])
rgm2s <- manifesto_reg_all_s(df_manifestos[
    term == "germany|teuton|prussia|fatherland|kaiser"])
rgm3s <- manifesto_reg_all_s(df_manifestos[
    term == "naval|navy|battleship|dreadnought|fleet"])
rgm4s <- manifesto_reg_all_s(df_manifestos[
    term == "national service league|navy league"])
table_manif_german_s <- stargazer(rgm1s[[1]], rgm1s[[2]], rgm2s[[1]], rgm2s[[2]],
    rgm3s[[1]], rgm3s[[2]], rgm4s[[1]], rgm4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, standardized.
        All models include constituency, party, and year fixed effects.
        ``German terms'' are ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,''
        and ``fatherland,'' ``Navy terms'' are ``navy,'' ``naval,''
        ``dreadnought,'' ``battleship,'' and ``fleet,''
        ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by county in parentheses.}",
    title = c("Effects of import competition on manifesto references to Germany
        and the naval race"),
    label = "table_manif_german_cl",
    notes.align = "l")
cat(table_manif_german_s, sep = '\n', file = "Tables/table_manif_german_cl.tex")

table_manif_german_s_full <- stargazer(rgm1s[[1]], rgm1s[[2]], rgm2s[[1]], rgm2s[[2]],
    rgm3s[[1]], rgm3s[[2]], rgm4s[[1]], rgm4s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.4\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, standardized.
        All models include constituency, party, and year fixed effects.
        ``German terms'' are ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,''
        and ``fatherland,'' ``Navy terms'' are ``navy,'' ``naval,''
        ``dreadnought,'' ``battleship,'' and ``fleet,''
        ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-19"),
    label = "table_manif_german_cl",
    notes.align = "l")
cat(table_manif_german_s_full, sep = '\n', file = "Tables with full output/table_a19.tex")



# A-41: effects of IPW on manifesto references to Germany and naval race, with
# exposure-robust SEs
rgm1 <- manifesto_reg_all(df_manifestos[term == "germany"])
rgm2 <- manifesto_reg_all(df_manifestos[
    term == "germany|teuton|prussia|fatherland|kaiser"])
rgm3 <- manifesto_reg_all(df_manifestos[
    term == "naval|navy|battleship|dreadnought|fleet"])
rgm4 <- manifesto_reg_all(df_manifestos[
    term == "national service league|navy league"])

models <- list(rgm1[[1]], rgm1[[2]], rgm2[[1]], rgm2[[2]],
    rgm3[[1]], rgm3[[2]], rgm4[[1]], rgm4[[2]])

f_stats <- sapply(models, get_f_stat)

table_manif_german <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``germany''", "German terms", "Navy terms", "Militarist groups"),
    column.separate = c(2, 2, 2, 2),
    add.lines = list(c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_manif_german_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022). Manifesto-level data aggregated to the industry level.
        Dependent variable is number of uses of
        specified term relative to total length of manifesto, standardized.
        All models include constituency, party, and year fixed effects.
        ``German terms'' are ``germany,'' ``kaiser,'' ``teuton,'' ``prussia,''
        and ``fatherland,'' ``Navy terms'' are ``navy,'' ``naval,''
        ``dreadnought,'' ``battleship,'' and ``fleet,''
        ``Militarist groups'' are ``national service league'' and
        ``navy league.'' Standard errors clustered by industry in parentheses.}",
    title = c("Effects of import competition on manifesto references to Germany
        and the naval race, exposure-robust standard errors"),
    label = "table_manif_german_ex",
    notes.align = "l")
cat(table_manif_german, sep = '\n', file = "Tables/table_manif_german_ex.tex")


# Table A-20 effects of IPW on references to immigration in Conservative campaign
# manifestos

ra_1s <- manifesto_reg_alien_s(df_manifestos[term == "immigra[un]t"])
ra_2s <- manifesto_reg_alien_s(df_manifestos[term == "alien"])
ra_3s <- manifesto_reg_alien_s(df_manifestos[term == "jew"])
ra_4s <- manifesto_reg_alien_s(df_manifestos[term == "foreigner"])
ra_5s <- manifesto_reg_alien_s(df_manifestos[
    term == "alien|immigra[un]t|jew|foreigner"])
table_manif_aliens_s <- stargazer(ra_1s[[1]], ra_1s[[2]], ra_2s[[1]], ra_2s[[2]],
    ra_3s[[1]], ra_3s[[2]], ra_4s[[1]], ra_4s[[2]], ra_5s[[1]], ra_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("as.factor", "const_frac_secondary", "const_frac_n_irish"),
    column.labels = c("``immigrant''", "``alien''", "``jew''", "``foreigner''",
        "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x", "x",
        "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Conservative
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by county in parentheses.}",
    title = c("Effect of local trade shocks on references to immigration in
        Conservative campaign manifestos"),
    label = "table_manif_aliens_cl",
    notes.align = "l")
cat(table_manif_aliens_s, sep = '\n',
    file = "Tables/table_manif_aliens_cl.tex")

table_manif_aliens_s_full <- stargazer(ra_1s[[1]], ra_1s[[2]], ra_2s[[1]], ra_2s[[2]],
    ra_3s[[1]], ra_3s[[2]], ra_4s[[1]], ra_4s[[2]], ra_5s[[1]], ra_5s[[2]],
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    #omit = c("as.factor", "const_frac_secondary", "const_frac_n_irish"),
    column.labels = c("``immigrant''", "``alien''", "``jew''", "``foreigner''",
        "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x", "x",
        "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x")),
    #keep.stat = c("n"),
    omit.stat = "ser",
    column.sep.width = "-10pt",
    #float.env = "sidewaystable",
    notes = "\\parbox[t]{0.4\\textwidth}{
        Manifesto-level regressions. Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Conservative
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by county in parentheses.}",
    title = c("Full regression output for Table A-20"),
    label = "table_manif_aliens_cl",
    notes.align = "l")
cat(table_manif_aliens_s_full, sep = '\n',
    file = "Tables with full output/table_a20.tex")


rm(ra_1s, ra_2s, ra_3s, ra_4s, ra_5s)

# Table A-42: effects of IPW on references to immigration in Conservative campaign
# manifestos, with exposure-robust SEs
ra_1 <- manifesto_reg_alien(df_manifestos[term == "immigra[un]t"])
ra_2 <- manifesto_reg_alien(df_manifestos[term == "alien"])
ra_3 <- manifesto_reg_alien(df_manifestos[term == "jew"])
ra_4 <- manifesto_reg_alien(df_manifestos[term == "foreigner"])
ra_5 <- manifesto_reg_alien(df_manifestos[
    term == "alien|immigra[un]t|jew|foreigner"])

models <- list(ra_1[[1]], ra_1[[2]], ra_2[[1]], ra_2[[2]],
    ra_3[[1]], ra_3[[2]], ra_4[[1]], ra_4[[2]], ra_5[[1]], ra_5[[2]])
f_stats <- sapply(models, get_f_stat)

table_manif_aliens <- stargazer(models,
    covariate.labels = c("$\\Delta\\textrm{IPW}_{1885}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    omit = c("Constant"),
    column.labels = c("``immigrant''", "``alien''", "``jew''", "``foreigner''",
        "All"),
    column.separate = c(2, 2, 2, 2, 2),
    add.lines = list(c("Initial immigrants x year", "x", "x", "x", "x", "x", "x",
        "x", "x", "x", "x"),
        c("Initial Mf x year", "", "x", "", "x", "", "x", "", "x",
        "", "x"),
        c("First stage F-stat", f_stats)),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    float.env = "sidewaystable",
    notes = "\\parbox[t]{0.6\\textwidth}{
        This table replicates the results of Table
        \\ref{table_manif_aliens_cl} using the aggregation
        and standard error calculation methods recommended by Borusyak, Hull,
        and Jaravel (2022).
        Manifesto-level data aggregated to the industry level.
        Dependent variable is number of uses of
        specified term relative to total length of manifesto, by Conservative
        candidates, standardized. All models include
        constituency and election fixed effects.
        Standard errors clustered by industry in parentheses.}",
    title = c("Effect of local trade shocks on references to immigration in
        Conservative campaign manifestos, exposure-robust standard errors"),
    label = "table_manif_aliens_ex",
    notes.align = "l")
cat(table_manif_aliens, sep = '\n',
    file = "Tables/table_manif_aliens_ex.tex")
rm(models)
rm(ra_1, ra_2, ra_3, ra_4, ra_5)

#%% Other tables
# Table 1: summary statistics
my_sums <- function(data, variable, new_name, category){
    #variable <- c(variable)

    var_vec <- unlist(data[, ..variable])
    var_vec <- na.omit(var_vec)
    var_n <- length(var_vec)
    var_sd <- sd(var_vec)
    var_mean <- mean(var_vec)
    var_min <- min(var_vec)
    var_max <- max(var_vec)
    return(c(category, new_name, var_n, var_mean, var_sd, var_min, var_max))
}
df_cross_section <- unique(df_s[!is.na(constituency_id), .(constituency_id, const_frac_secondary,
    const_frac_n_irish)])
sum_const_frac_secondary <- my_sums(df_cross_section, "const_frac_secondary",
    "Manufacturing share 1881", "Constituency")
sum_const_frac_n_irish <- my_sums(df_cross_section, "const_frac_n_irish",
    "Immigrant share 1881", "Constituency")
sum_shock_diff_w <- my_sums(df_s[unopposed == 0], "shock_diff_w", "$\\Delta\\textrm{IPW}_{1885}$",
    "Constituency x election")
sum_frac_secondary <- my_sums(df_s[year %in% c(1880, 1890, 1900, 1910)],
    "frac_secondary", "Manufacturing share",
    "Constituency x census year")
sum_frac_unskilled <- my_sums(df_s[year %in% c(1880, 1890, 1900, 1910)],
    "frac_unskilled", "Unskilled jobs share", "Constituency x census year")
sum_frac_vagrant <- my_sums(df_s[year %in% c(1880, 1890, 1900, 1910)],
    "frac_vagrant", "Vagrant share", "Constituency x census year")
sum_hiscam <- my_sums(df_s[year %in% c(1880, 1890, 1900, 1910)],
    "hiscam", "Average economic status", "Constituency x census year")
sum_cons <- my_sums(df_s[!is.na(shock_diff_w) & unopposed == 0], "cons_share", "Conservative vote share",
    "Constituency x election")
sum_lib <- my_sums(df_s[!is.na(shock_diff_w) & unopposed == 0], "lib_share", "Liberal vote share",
    "Constituency x election")
sum_lab <- my_sums(df_s[!is.na(shock_diff_w) & unopposed == 0], "lab_share", "Labour vote share",
    "Constituency x election")
sum_frac_category <- my_sums(unique(df[year == 1900,
    .(constituency_id, category_uk, const_frac_category)]),
    "const_frac_category", "Industry share", "Constituency x industry")
sum_frac_category
44080 - 43985
43985 / 95
df_s[!(constituency_id %in% df$constituency_id) & !is.na(const_frac_secondary),
    .(constituency_id, const_frac_secondary)]


sum_diff_percap <- my_sums(unique(df[,
    .(category_uk, year, diff_percap_w)]), "diff_percap_w",
        "$\\Delta\\textrm{IPW}_{1885}$", "Industry x year")
sum_fd_hiscam <- my_sums(df_s[year %in% c(1890, 1900, 1910)],
    "fd_hiscam", "$\\Delta$Average economic status",
        "First difference constituency x year")
sum_fd_ln_frac_unskilled <- my_sums(df_s[year %in% c(1890, 1900, 1910)],
    "fd_ln_frac_unskilled", "$\\Delta$ln unskilled jobs share",
        "First difference constituency x year")
sum_fd_ln_frac_vagrant <- my_sums(df_s[year %in% c(1890, 1900, 1910)],
    "fd_ln_frac_vagrant", "$\\Delta$ln vagrant share",
        "First difference constituency x year")
sum_fd_shock <- my_sums(df_s[year %in% c(1890, 1900, 1910)],
    "fd_shock_census_w", "$\\Delta\\textrm{IPW}_t$ ",
        "First difference constituency x year")
df_sum_table <- rbind.data.frame(
    sum_const_frac_secondary,
    sum_const_frac_n_irish,
    sum_frac_category,
    sum_frac_secondary,
    sum_frac_vagrant,
    sum_frac_unskilled,
    sum_hiscam,
    sum_fd_shock,
    sum_fd_ln_frac_vagrant,
    sum_fd_ln_frac_unskilled,
    sum_fd_hiscam,
    sum_shock_diff_w,
    sum_cons,
    sum_lib,
    sum_lab
    )
df_sum_table[, 3:7] <- apply(df_sum_table[, 3:7], 2, as.numeric)
colnames(df_sum_table) <- c("Category", "Variable", "N", "Mean", "SD", "Min", "Max")


df_sum_table

print(df_sum_table)
sum_table <- kbl(df_sum_table[, 2:7], "latex", booktabs = T, digits = 3, format.args = list(big.mark = ","),
    caption = "Summary statistics", escape = FALSE) %>%
    kable_styling(position = "center") %>%
    pack_rows(index = c("Constituency" = 2,
        "Constituency x Industry" = 1,
        "Constituency x census year" = 4,
        "First difference constituency x census year" = 4,
        "Constituency x election year" = 4),
        bold = F, italic = T)
sum_table
cat(sum_table, sep = "/n", file = "Tables/sum_stats.tex")


# A-1: employment and 1885-1910 shock by industry
df_categories <- df[year == 1910, .(
    emp_1881 = sum(const_frac_category * const_pop)),
    by = .(category_uk, diff_percap_w)]
setorder(df_categories, -emp_1881)

df_categories <- df_categories[category_uk != "others",
    .(Industry = str_to_title(category_uk),
    `1881 Employment` = emp_1881,
    `$\\Delta$IPW (1885-1910)` = diff_percap_w)]
df_categories
category_table <- kbl(df_categories, "latex", booktabs = T, digits = 3,
    format.args = list(big.mark = ","),
    caption = "Industry categories", escape = FALSE,
    longtable = TRUE) %>%
    kable_styling(position = "center", latex_options = "hold_position")
cat(category_table, sep = "/n", file = "Tables/category_table.tex")


# A-44: regression version of Fig 5
rmf1 <- felm(cons_share ~ treat_1910 : as.factor(year) | year | 0 | county_name ,
    df_s)
summary(rmf1)
rmf2 <- felm(cons_share ~ treat_1910 : as.factor(year) | year | 0 | county_name,
    df_match)
summary(rmf2)
table_match_plot <- stargazer(rmf1, rmf2,
    covariate.labels = c(
        "$\\Delta\\textrm{IPW}_{1900}$, 1900--1910 above median \\\\$\\times \\text{election} = 1885$",
        "$\\times \\text{election} = 1886$", "$\\times \\text{election} = 1892$",
        "$\\times \\text{election} = 1895$", "$\\times \\text{election} = 1900$",
        "$\\times \\text{election} = 1906$", "$\\times \\text{election} = 1910\\text{J}$",
        "$\\times \\text{election} = 1910\\text{D}$"),
    dep.var.labels.include = F,
    dep.var.caption = "",
    add.lines = list(c("Matched panel", "", "x")),
    keep.stat = c("n"),
    #omit.stat = "ser",
    column.sep.width = "-10pt",
    notes = "\\parbox[t]{0.5\\textwidth}{
        This table reports the results of regressions of Conservative vote share
        on a dummy variable which takes a value of 1 if the change in imports
        per worker between 1900 and 1910 was above the median, and zero if
        otherwise, interacted with election fixed effects. This provides a
        regression-based complement to Figure 5. Model (1) is estimated using the
        full panel, (2) is estimated using a panel matched on 1885, 1892, and 1900
        Conservative vote shares.
        All models include election fixed effects.
        Standard errors clustered by county in parentheses.}",
    title = c(
        "Evolution of Conservative voting by 1900--1910 change in imports, with
        full and matched panel"),
    label = "table_match_plot",
    notes.align = "l"
)
cat(table_match_plot, sep = '\n',
    file = "Tables/table_match_plot.tex")


#%% Plots
# figure 1: German imports over time

df_f1 <- fread("Data/_Data for figures/figure_1.csv")
ggplot(df_f1, aes(x = year, y = imports_combined)) +
    geom_line() +
    labs(x = element_blank(),
        y = "Imports from Germany, the Netherlands\nand Belgium (Million pounds)",
        title = "German imports by year")
ggsave("Figures/figure_1.png", width = 5, height = 3,
    dpi = 300)

# figure 2
df_import_cat <- fread("Data/_Data for figures/figure_2.csv")
head(df_import_cat)

stone_and_ore <- c("copper ore", "iron ore", "lead ore", "silver ore", "slate",
    "stone", "tin ore", "gold ore", "sand flint clay gravel chalk")
sheet_metal <- c("sheet gold silver", "sheet lead", "sheet zinc", "sheet copper",
    "sheet iron and steel", "sheet other metals", "tin")
metal_goods <- c("brass manufactures", "buttons", "implements and tools",
    "machinery", "copper manufactures", "hardware and cutlery", "iron manufactures",
    "zinc manufactures")
food_drink_tobacco <- c("beer", "dairy", "fish", "meat", "mustard vinegar spice pickle",
    "refined sugar", "spirits", "tobacco manufactures", "chocolate",
    "soft drinks", "jams and sweets")
textile_yarn <- c("wool yarn", "cotton yarn", "silk yarn")
textiles_apparel <- c("hats", "lace", "apparel and haberdashery",
    "cordage",  "gloves", "jute manufactures", "linen",
    "mats", "plaiting of straw",  "wool manufactures",
    "cotton manufactures", "embroidery", "silk manufactures")
chemical_products <- c("alkali", "asbestos", "chemicals", "dyes and paints",
    "manure", "oil seed and oil cake", "candles and grease",
    "floor cloth and oil cloth", "glue", "matches", "soap")
skins_and_hides <- c("skins and furs", "leather")
paper_goods <- c("paper", "printed matter")
manufactured_goods <- c("arms and ammunition", "clocks and watches", "electricals",
    "umbrellas and sticks", "bicycles", "carriages", "motor cars",
    "scientific instruments", "toys", "lamps", "machinery", "musical instruments")
leather_goods <- c("gloves", "shoes", "leather manufactures")
rubber_goods <- c("rubber", "gum", "waterproof goods")
other_goods <- c("art", "artificial flowers", "bristles and brushes", "cement",
    "fancy goods", "feathers", "hay", "jewellery",
    "tobacco pipes")
df_import_cat[category_uk %in% c(stone_and_ore, "coal coke and patent fuel"),
    big_category := "Coal, stone and ore"]
df_import_cat[category_uk %in% sheet_metal, big_category := "Sheet metal"]

df_import_cat[category_uk %in% food_drink_tobacco,
    big_category := "Food, drink and tobacco"]
df_import_cat[category_uk %in% textiles_apparel,
    big_category := "Textile products and apparel"]
df_import_cat[category_uk %in% textile_yarn,
    big_category := "Textile yarn"]
df_import_cat[category_uk %in% c(chemical_products, rubber_goods),
    big_category := "Chemicals and related products"]

df_import_cat[category_uk %in% c(skins_and_hides, leather_goods),
    big_category := "Hides and leather"]
df_import_cat[category_uk %in% paper_goods,
    big_category := "Paper products"]
df_import_cat[category_uk %in% c(manufactured_goods, metal_goods), big_category := "Manufactured goods"]


df_import_cat[category_uk %in% c(paper_goods, "glass", "wood products"),
    big_category := "Paper, glass, and wood"]

df_import_cat[category_uk %in% c(other_goods, "china and earthenware"),
    big_category := "Other products"]
df_import_big_cat <- df_import_cat[, .(imports = sum(value_tot / 1000000)),
    by = .(big_category, year)]

ggplot(df_import_big_cat[],
    aes(x = year, y = imports, color = big_category)) +
    geom_line() +
    geom_text_repel(data = df_import_big_cat[year == 1910],
        aes(x = Inf, y = imports, label = big_category), color = "black") +
    labs(x = element_blank(), color = element_blank()) +
    scale_color_brewer(palette = "RdYlBu") +
    #scale_color_viridis_d(option = "plasma") +
    theme(legend.position = "none") +
    labs(x = element_blank(),
        y = "Imports from Germany, the Netherlands and Belgium (million pounds)",
        title = "German imports by category, 1880-1910")
ggsave("Figures/figure_2.png", width = 10, height = 6,
    dpi = 300)

# figure 3
df_f3 <- fread("Data/_Data for figures/figure_3.csv")

ggplot(df_f3, aes(x = year, y = total_welfare)) +
    geom_line() +
    labs(x = element_blank(), y = "Social welfare spending (million pounds)",
        title = "Social welfare spending")
ggsave("Figures/figure_3.png", width = 5, height = 3, dpi = 300)



# Figure 5
# note -- must run the section for Table 3 first, which generates the matched
# dataset
df_s[treat_1910 == 0, treat_1910_text := "below median"]
df_s[treat_1910 == 1, treat_1910_text := "above median"]
df_match[treat_1910 == 0, treat_1910_text := "below median"]
df_match[treat_1910 == 1, treat_1910_text := "above median"]
plot_vote1 <- ggplot(df_s[!is.na(treat_1910) & year >= 1885],
  aes(x = year, y = cons_share, color = treat_1910_text)) +
  geom_jitter(alpha = 0.2) +
  geom_smooth(method = "loess") +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Conservative vote share by 1910 trade shock",
       x = element_blank(),
       y = "Conservative vote share",
       color = expression(Delta*"IPW, 1900-1910"))
plot_vote2 <- ggplot(df_match[!is.na(treat_1910) & year >= 1885],
  aes(x = year, y = cons_share, color = treat_1910_text,
  weight = weight)) +
  geom_jitter(alpha = 0.2) +
  ylim(c(0, 1)) +
  geom_smooth(method = "loess") +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Conservative vote share, matched panel",
       x = element_blank(),
       y = "Conservative vote share",
       color = expression(Delta*"IPW, 1900-1910"))

plots_vote <- grid.arrange(plot_vote1, plot_vote2, ncol = 2)


ggsave("Figures/figure_5.png", plots_vote, width = 10, height = 6,
    dpi = 300)



# figure 6
df_f6 <- fread("Data/_Data for figures/figure_6.csv")
head(df_f6)
df_f6[, term := str_replace_all(term, "\"\"", "\"")]
ggplot(df_f6[between(year, 1880, 1910)], aes(x = year, y = references, color = term)) +
    geom_line() +
    geom_text_repel(data = df_f6[year == 1910],
        aes(x = Inf, y = references, label = term), color = "black", nudge_y = c(0, 25, 40)) +
    labs(x = element_blank(), color = element_blank()) +
    scale_color_brewer(palette = "Set1") +
    theme(legend.position = "none") +
    labs(x = element_blank(),
        y = "References in the Times",
        title = "References to unemployment in the Times newspaper,\n1880-1910")
ggsave("Figures/figure_6.png", width = 5, height = 3, dpi = 300)


# figure A-1
df_fa1 <- fread("Data/_Data for figures/figure_a1.csv")
head(df_fa1)
df_fa1[, type := factor(type,levels = c("Port", "Origin"))]
df_fa1[, country := factor(country,levels = c("Germany", "Belgium", "Netherlands"))]

ggplot(df_fa1, aes(x = type, y = imports, fill = country)) +
    geom_bar(position = "dodge", stat = "identity") +
    labs(x = "Method of classification", y = "Imports (million pounds)",
        title = "Comparing import country attribution systems for 1910",
        fill = "Country") +
    scale_fill_brewer(palette = "Set1")
ggsave("Figures/figure_a1.png", width = 7.5, height = 5,
    dpi = 300)

# figure A-2
df_trad_hist <- fread("Data/_Data for figures/figure_a2.csv")

# leading destinations of German exports
df_trad_hist[iso_o %in% c("DEU", "BEL", "NLD", "PRUS") & year == 1900, .(FLOW = sum(FLOW)),
    by = .(iso_d, year)][order(-FLOW)]

df_trad_hist_plot <- df_trad_hist[iso_o %in% c("DEU", "BEL", "NLD") & between(year, 1880, 1910)
        & iso_d %in% c("GBR", "FRA", "AUTHUN", "USA", "USSR", "CHE", "SWE", "DNK", "ITA"),
        .(FLOW = sum(FLOW, na.rm = T)),
        by = .(iso_d, year)]
df_trad_hist_plot[iso_d == "USSR", iso_d := "RUS"]



plot_german_export <- ggplot(df_trad_hist_plot, aes(x = year, y = FLOW / 1e6, group = paste(iso_d), color = iso_d == "GBR")) +
    geom_line() +
    geom_text_repel(data = df_trad_hist_plot[year == 1910],
        aes(x = Inf, y = FLOW / 1e6, label = paste(iso_d))) +
    theme(legend.position = "none") +
    scale_color_manual(values = c("#000000", brewer.pal(n = 9, name = "Set1")[1])) +
    labs(x = element_blank(), color = element_blank(),
        y = "Exports from Germany, Netherlands and Belgium (million pounds)",
        title = "German exports, 1880-1910")
df_trad_hist[year == 1900 & iso_d == "GBR"][order(-FLOW)][1:20]
df_trad_hist_plot2 <- df_trad_hist[iso_d %in% c("GBR") & between(year, 1880, 1910)
        & iso_o %in% c("USA", "DEU", "BEL", "NLD", "FRA", "GBRIND", "AUS",
            "CAN", "USSR", "ESP", "DNK"),
        .(FLOW = sum(FLOW, na.rm = T)),
        by = .(iso_o, year)]
df_trad_hist_plot2[iso_o %in% c("NLD", "BEL", "DEU"), iso_o := "DEU"]
df_trad_hist_plot2[iso_o == "USSR", iso_o := "RUS"]
df_trad_hist_plot2[iso_o == "GBRIND", iso_o := "IND"]
df_trad_hist_plot2 <- df_trad_hist_plot2[, .(FLOW = sum(FLOW, na.rm = T)),
    by = .(iso_o, year)]
plot_uk_import <- ggplot(df_trad_hist_plot2, aes(x = year, y = FLOW / 1e6, group = paste(iso_o), color = iso_o == "DEU")) +
    geom_line() +
    geom_text_repel(data = df_trad_hist_plot2[year == 1910],
        aes(x = Inf, y = FLOW / 1e6, label = paste(iso_o))) +
    theme(legend.position = "none") +
    scale_color_manual(values = c("#000000", brewer.pal(n = 9, name = "Set1")[1])) +
    labs(x = element_blank(), color = element_blank(),
        y = "UK imports (million pounds)",
        title = "UK imports, 1880-1910")
trade_grid <- grid.arrange(plot_german_export, plot_uk_import, ncol = 2)

ggsave("Figures/figure_a2.png", trade_grid, width = 10, height = 6,
    dpi = 300)

# figure 4
# Note: this figure requires use of Great Britain Historical GIS shapefiles
# for English and Welsh constituencies, 1885-1918, which must be obtained from
# www.VisionofBritain.org.uk

# replace these file locations with the location of your version of the
# shapefiles
england_shapefiles_string <- paste0("~/Dropbox (Personal)/Trade shocks 1900s/",
    "Data/Constituency Boundaries/Boundaries 1885/",
    "E1885_constituencies/E1885_1918_Constituencies.shp")
wales_shapefiles_string <- paste0("~/Dropbox (Personal)/Trade shocks 1900s/",
    "Data/Constituency Boundaries/Boundaries 1885/",
    "W1885_constituencies/W1885_1918_Constituencies.shp")

sf_const_e <- st_read(england_shapefiles_string)

sf_const_w <- st_read(wales_shapefiles_string)
head(sf_const_w)

sf_const <- rbind(sf_const_e, sf_const_w)

df_s_map <- merge(sf_const, df_s, by = "g_unit")
print(colnames(df))




bin_quintiles <- function(vec){
    quintiles <- quantile(vec, c(0, 0.2, 0.4, 0.6, 0.8, 1))
    binned_vec <- cut(vec, quintiles, labels = F)
    return(binned_vec)
}

df_s_map$shock_diff_w_bin[df_s_map$year == 1910] <- bin_quintiles(
    df_s_map$shock_diff_w[df_s_map$year == 1910])


map_shock <- ggplot(subset(df_s_map, year == 1910 & !is.na(shock_diff_w_bin)),
    aes(fill = as.factor(shock_diff_w_bin))) +
    geom_sf(size = 0) +
    scale_fill_brewer() +
    labs(fill = expression(Delta*"IPW quintile"),
        title = "Geographic distribution of the trade shock, 1910\n")

r_shock_resid <- felm(shock_diff_w ~ const_frac_secondary | 0,
    data = subset(df_s_map, year == 1910))
df_s_map$shock_resid[df_s_map$year == 1910] <- r_shock_resid$residuals
df_s_map$shock_resid_bin[df_s_map$year == 1910] <- bin_quintiles(
    df_s_map$shock_resid[df_s_map$year == 1910]
)

map_shock_resid <- ggplot(subset(df_s_map, year == 1910 & !is.na(shock_resid_bin)),
    aes(fill = as.factor(shock_resid_bin))) +
    geom_sf(size = 0) +
    scale_fill_brewer() +
    labs(fill = expression(Delta*"IPW residualized quintile"),
        title = paste0("Geographic distribution of the trade shock, 1910\n",
            "controlling for initial manufacturing"))

map_shock_resid

maps <- grid.arrange(map_shock, map_shock_resid, ncol = 2)


ggsave("Figures/figure_4.png", maps, width = 10, height = 6,
    dpi = 300)
